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Summary 


There  is  growing  interest  in  validating  existing  numerical  methods  and  physical  models 
employed  in  currently  state-of-the-art  computational  fluid  dynamics  (CFD)  codes.  The 
numerical  methods  and  physical  models  that  are  under  examination  and  are  the  subject  of 
validation  studies  pertain  to  the  simulation  of  hypersonic  low  Reynolds  number  flows  that 
are  in  the  continuum  and  the  near  continuum  limit.  This  work  was  based  on  the  framework 
that  we  established  for  the  validation  of  numerical  methods  for  inert  hypersonic  flows  under 
Air  Force  support  in  recent  years. 

In  the  past,  we  simulated  experiments  of  hypersonic  laminar  double-cone  flows  that 
were  performed  at  the  Large  Energy  National  Shock  (LENS)  facility.  The  experiments 
were  conducted  as  part  of  a  validation  effort  for  continuum  and  par  tide- based  codes.  The 
double-cone  flow  was  chosen  because  it  exhibits  strong  viscous/inviscid  and  shock  interac¬ 
tions,  and  we  have  shown  that  the  separation  zone  that  forms  at  the  cone-cone  juncture 
is  sensitive  to  both  the  amount  of  numerical  dissipation  in  the  numerical  scheme  and  the 
amount  of  energy  that  participates  in  the  dissociation  process  and  other  chemical  relax¬ 
ation  phenomena  (real-gas  effects).  The  double-cone  flows  are  very  challenging  to  compute 
even  with  present  day  computing  power.  As  a  result  of  the  previous  effort,  a  large  number 
of  experiments  were  performed  in  air  under  high  enthalpy  conditions.  The  experiments 
were  performed  in  different  facilities  at  LENS,  and  data  exists  for  a  range  of  frce-stream 
Mach  and  Reynolds  numbers  at  various  level  of  free-stream  enthalpy.  Simulations  of  these 
experiments  that  employ  standard  chemical  kinetics  models  showed  poor  agreement  with 
the  experimental  data.  Agreement  worsens  as  the  free-stream  enthalpy  is  increased.  Most 
of  present  work  focuses  on  the  analysis  of  the  experimental  conditions.  Effort  was  put 
forth  to  analyze  these  high  enthalpy  experiments  by  employing  the  most  up-to-date  phys¬ 
ical  models  of  high  enthalpy  air  chemistry. 

In  particular,  we  used  a  vibrational  state-specific  model  for  vibrational  relaxation  in 
the  LENS  nozzle  flow.  We  also  considered  a  wide  variety  of  thermo-chemical  models  and 
chemical  kinetics  models  for  rapid  expansion  of  air  in  a  hypersonic  nozzle.  To  date,  these 
modeling  approaches  have  been  unsuccessful  for  high  enthalpy  air.  We  find  that  once  the 
equilibrium  composition  of  the  reservoir  gas  has  an  appreciable  level  of  oxygen  atoms,  we 
obtain  poor  agreement  with  measurements  in  the  LENS  facility.  This  occurs  above  about 
5  MJ/kg  in  air.  We  are  continuing  to  work  on  this  problem,  and  this  report  documents 
what  we  have  tried  to  date. 


I.  Introduction 


1.  Motivations 

There  has  been  a  great  deal  of  interest  in  developing  prototype  re-entry  vehicles  and 
hypersonic  air-breathing  vehicles  for  access  to  space.  These  vehicles  will  be  designed  to 
operate  in  the  high  enthalpy  regime,  where  the  (low  exhibits  high  enthalpy  or  “real-gas”  ef¬ 
fects.  Real -gas  effects  are  triggered  by  the  high  temperatures  encountered  in  high  enthalpy 
flows.  These  effects  include  excitation  of  the  internal  energy  modes  of  the  gas,  as  well  as 
dissociation  of  polyatomic  species  and  ionization.  These  physical  phenomena  occur  over  a 
range  of  time  scales,  and  the  flow  is  frequently  in  chemical  and  thermal  nonequilibrium. 
Additionally,  depending  on  the  shape  of  the  vehicle  and  the  flight  conditions,  shock/shock 
interactions  and  viscous-inviscid  interactions  will  likely  occur  during  flight  at  hypersonic 
speeds.  The  combined  effects  of  these  phenomena  and  real- gas  chemistry  may  dramatically 
alter  the  flight  characteristics  of  the  vehicle.  Therefore,  all  of  these  effects  must  be  taken 
into  consideration  in  the  development  stage  of  future  hypersonic  vehicles. 

The  development  and  design  of  such  vehicles  is  aided  by  the  use  of  experimentation 
and  numerical  simulation.  There  are  difficulties  associated  with  both  approaches.  It  is 
difficult  to  produce  ground-based  hypersonic  flow  experiments  that  are  realistic  in  terms 
of  the  enthalpy  and  composition  of  the  gas.  It  is  also  difficult  to  perform  accurate  numerical 
simulations,  as  the  physics  involved  is  not  fully  understood. 

2.  Double-Cone  Flows 

Flows  generated  by  double-cone  geometries  placed  in  hypersonic  free-streams  have 
been  studied  in  the  past.  Earlier  studies  of  these  flows  focused  primarily  on  testing  and 
improving  nonequilibrium  chemistry  models.  This  is  because  the  double-cone  flows  are 
sensitive  to  high  enthalpy  effects.  These  studies  were  largely  unsuccessful.  In  more  recent 
work,  we  studied  double-cone  flows  for  the  purpose  of  validating  numerical  methods.  This 
was  the  first  step  in  a  “building- block”  approach  to  assess  the  ability  of  computational 
fluid  dynamics  codes  to  reproduce  complex  hypersonic  flows  from  the  low  enthalpy  to  the 
high  enthalpy  regime.  The  effort  was  proven  successful  in  terms  of  validating  the  numerical 
methods  by  making  comparisons  with  experimental  measurements  of  double-cone  flows  in 
inert  environments  at  low  enthalpy  conditions. 

The  flow  generated  by  a  double-cone  configuration  placed  in  a  hypersonic  free-stream, 
exhibits  strong  viscous-inviscid  and  shock/shock  interactions.  The  oblique  shock  generated 
by  the  first,  shallow-angle  cone,  interacts  with  the  shock  generated  by  the  second  cone. 
When  the  second  cone  angle  is  larger  than  a  critical  angle,  this  shock  is  detached,  and  the 
interaction  is  very  strong.  A  schematic  of  the  flow  generated  by  a  double-cone  in  a  nitrogen 
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free-stream  at  a  Mach  number  of  about  12  is  shown  in  Fig,  1.  In  this  figure,  the  model 
lias  a  first  cone  half- angle  of  25°  and  a  second  cone  half- angle  of  55°.  The  shock  structure 
and  subsonic  regions  are  outlined,  and  the  shear  layers  are  visualized.  The  oblique  shock 
interacts  with  the  bow  shock  generated  by  the  second  cone,  and  a  transmitted  shock  is 
formed,  which  impinges  on  the  surface  downstream  of  the  cone-cone  juncture.  The  adverse 
pressure  gradient  at  the  cone-cone  juncture  causes  the  boundary  layer  to  separate.  The 
separated  region  that  forms  at  the  corner  creates  its  own  shock,  which  interacts  with  the 
bow  shock.  This  shifts  the  point  of  interaction,  which  in  turn  alters  the  separation  zone. 
This  process  feeds  back  on  itself  until  the  flow  reaches  steady  state  -  if  a  steady  state 
exists.  A  supersonic  jet  is  formed  along  the  surface  of  the  second  cone  downstream  of  the 
impingement  point,  and  it  undergoes  a  series  of  isentropic  compressions  and  expansions. 
The  flow  behind  the  strong  bow  shock  is  subsonic,  and  therefore  the  shape  of  the  jet 
influences  the  shock  structure, 


Figure  1.  Schematic  diagram  of  a  typical  flow- field  generated  by  a  25°-55°  double¬ 
cone  geometry  placed  in  a  Mach  12  free-stream.  The  flow  is  from  left  to  right,  and 
a  slice  is  shown  of  flow  that  is  circularly  symmetric  about  the  axis  of  the  geometry. 
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Under  high  enthalpy  conditions,  the  chemical  processes  that  take  place  in  the  How 
are  tightly  coupled  with  the  fluid  motion  through  energy  exchange  mechanisms.  In  high 
enthalpy  double-cone  flows,  dissociation  takes  place  both  in  the  high  temperature  region 
behind  the  bow  shock,  and  inside  the  separation  zone.  In  these  regions,  the  highly  ener¬ 
getic  particles  undergo  collisions  and  lose  some  of  their  energy  by  causing  dissociation  of 
polyatomic  species.  By  this  process,  energy  is  removed  from  the  bulk  flow  and  stored  in 
the  form  of  chemical  energy.  The  rate  at  which  this  process  takes  place  greatly  impacts 
the  mean  flow.  Thus,  the  size  of  the  embedded  region  of  separation  is  sensitive  to  chem¬ 
ical  effects  in  the  flow.  Additionally,  the  size  of  the  separation  zone  may  be  measured 
experimentally  by  surface  measurements  such  as  surface  pressure  and  heat  transfer  rate. 
Thus,  we  can  make  direct  comparisons  between  numerical  predictions  and  experimental 
measurements, 

3.  Summary  of  Previous  Work 

We  have  studied  extensively  hypersonic  double-cone  flows  with  and  in  the  absence 
of  flowfield  chemistry.  In  earlier  work  related  to  CFD  code  validation,  we  analyzed  ex¬ 
periments  performed  at  the  Large  Energy  National  Shock  (LENS)  facility.  We  studied 
low  enthalpy  laminar  double-cone  flows  and  showed  that  because  these  experiments  were 
performed  at  low  pressure,  there  was  vibrational  energy  frozen  in  the  free-stream,  which  in 
turn  reduced  the  amount  of  kinetic  energy  that  was  used  to  infer  the  free-stream  conditions 
via  a  Pitot  pressure  measurement.  This  resulted  in  erroneously  inferring  the  free-stream 
conditions,  which  resulted  in  seeing  discrepancies  between  the  experimental  measurements 
and  the  numerical  predictions.  When  we  accounted  for  that  effect  in  numerical  simulations, 
and  also  augmented  the  boundary  conditions  formulation  to  account  for  the  slip  effects  due 
to  the  low  pressures  of  the  experiment,  we  observed  excellent  agreement  between  calcu¬ 
lations  and  experiments  for  a  number  of  cases.  Fig.  2  plots  the  computed  and  measured 
heat  transfer  rates  to  the  double-cone  model  for  low  enthalpy  nitrogen  flows  at  Mach  9 
(Run  28)  and  Mach  12  (Run  35).  These  calculations  were  performed  at  the  newly  inferred 
free-stream  conditions  and  with  the  slip  model  for  the  internal  energy  modes  at  the  wall. 
Fig.  3  plots  the  computed  and  measured  heat  transfer  rates  to  the  double-cone  model  for 
low  enthalpy  Mach  15  nitrogen  flows  at  very  low  pressures.  Run  4  is  referred  to  as  “1/2 
density  case"  and  Run  5  as  “1/3  density  case"  and  these  designations  are  with  reference 
to  Run  35.  These  two  experimental  conditions  were  designed  to  provide  a  test-case  for 
the  more  computationally  intensive  DSMC  methods.  Agreement  is  excellent  between  the 
CFD  predictions  and  the  experimental  measurements  for  these  cases. 

We  observed  good  agreement  for  flows  of  nitrogen  at  moderate  enthalpy  conditions 
(free-stream  specific  enthalpy  of  5  MJ/kg),  and  the  results  are  plotted  in  Fig.  4,  which 
plots  the  predicted  heat  transfer  to  the  wall  along  witli  experimental  data  for  Mach  10 
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How.  The  result  obtained  under  the  computed  conditions  is  very  similar  to  the  original 
predicted  result.  The  computed  conditions  were  obtained  by  performing  calculations  of 
the  nozzle  flow  with  CFD  at  the  conditions  of  the  experiment.  The  computed  conditions 
were  extracted  from  the  CFD  simulation  at  the  location  where  the  model  is  situated. 


Figure  2.  Computed  heat  transfer  rate  for  low  enthalpy  nitrogen  at  Mach  9  (Run 
28)  and  Mach  12  (Run  35). 


Figure  3.  Computed  heat  transfer  rate  to  the  wall  for  the  nominally  Mach  15 
nitrogen  flows  Run  4  (“1/2  density  case”)  and  Run  5  (“1/3  density  case”). 
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Figure  4,  Computed  heat  transfer  rate  to  the  wall  for  Mach  10  nitrogen  flow  at 
moderate  enthalpy  conditions  (Run  42). 


Figure  5.  Computed  heat  transfer  rate  to  the  wall  for  Mach  10  nitrogen  flow  at 
moderate  enthalpy  conditions  and  higher  free-stream  Reynolds  number  (Run  46). 
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Figure  6.  Computed  heat  transfer  rate  to  the  wall  for  Mach  10  air  flow  at  10 
MJ/kg  free-stream  specific  enthalpy. 

Calculations  of  double-cone  flows  in  high  enthalpy  air  showed  poor  agreement.  Fig.  C 
shows  the  computed  heat  transfer  rate  results  for  three  simulations  of  Mach  10  air  flow 
at  10  MJ/kg  free-stream  specific  enthalpy.  The  calculation  under  nominal  conditions  does 
not  predict  the  separation  zone  size,  and  the  heat  flux  to  the  wall  is  underpredicted  even 
in  the  attached  region  of  the  flow  along  the  first  cone.  The  peak  in  heat  transfer  is  not 
predicted.  Bounding  calculations  using  the  fictitious  super-catalytic  wall  condition  showed 
only  a  small  difference  in  heat  flux.  We  would  assume  that  any  agreement  with  the  data 
in  this  case  would  be  spurious.  Additional  calculations  of  the  double-cone  flow  under 
conditions  of  forced  chemical  equilibrium  (very  fast  reaction  rates  in  the  computed  nozzle 
flow)  and  forced  vibrational  equilibrium  do  not  show  improvement. 

Many  more  air  cases  were  simulated  and  similar  discrepancies  were  observed.  Figure 
7  plots  the  heat  transfer  rate  to  the  wall  for  three  air  cases  (Run  59  at  4.5  MJ/kg,  Run 
54  at  10.4  MJ/kg,  and  Run  G5  at  15.2  MJ/kg)  which  are  arranged  in  ascending  order  of 
free-stream  enthalpy.  We  see  that  as  the  enthalpy  increases  agreement  worsens. 
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Figure  7.  Predicted  and  measured  heat  transfer  rate  for  Run  59,  Run  54  and  Run  65 
(high  enthalpy  air  cases  ordered  in  terms  of  free-stream  enthalpy). 
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4.  Scope  of  the  Present  Work 

In  this  work,  we  use  computational  fluid  dynamics  to  study  hypersonic  double-cone 
experiments  performed  in  the  LENS  facility  at  high  enthalpy  air.  Because  wc  are  unable 
to  explain  and  correct  for  the  discrepancies  that  we  observe  between  numerical  predictions 
and  the  experimental  measurements,  we  study  the  experiments  using  advanced  chemistry 
models.  We  show  results  for  high  enthalpy  air  double-cone  flows  for  range  of  conditions. 
These  simulations  utilize  standard  chemical  kinetics  models  that  are  being  used  presently 
in  hypersonic  CFD  applications.  We  also  make  use  of  results  from  recent  work  on  non- 
intrusive  diagnostics  of  the  flow  inside  the  nozzle.  We  compare  the  predicted  nitric  oxide 
(NO)  content  at  the  nozzle  exit  plane  with  what  is  inferred  from  laser  absorption  measure¬ 
ments  made  at  LENS. 

The  following  two  sections  of  the  report  are  included  for  completeness.  Section  2 
concerns  the  governing  equations  that  are  solved  by  the  computational  fluid  dynamics 
code.  The  standard  models  that  are  employed  in  CFD  are  outlined.  Augmentations  to 
the  baseline  models  to  accommodate  for  more  advanced  flow  physics  are  discussed  in  the 
relevant  sections.  We  present  the  governing  equations  for  the  continuum  flow  of  a  gas  in 
vibrational  and  chemical  nonequilibrium.  The  Navier-Stokes  equations  for  a  chemically 
reacting  and  vibrationally  relaxing  air  mixture  are  presented,.  The  numerical  treatment  of 
the  equations  is  presented  in  section  3.  We  use  standard  numerical  techniques  which  are 
popular  among  researchers  and  have  been  widely  used  for  practical  applications.  It  should 
be  noted  that  we  have  validated  the  use  of  these  models  for  the  double-cone  flow  in  earlier 
work.  We  are  confident  that  with  the  numerical  treatment  of  choice  and  with  our  previous 
experience  in  grid  generation  for  these  problems  we  are  adequately  resolving  the  fluid  flow 
physics. 

Section  four  explores  the  possibility  of  employing  a  detailed  vibrational  relaxation  and 
diatomic  molecule  dissociation  model  based  on  a  recent  state-specific  modeling  formula¬ 
tion.  The  work  described  analyzes  the  sensitivity  of  the  free-stream  conditions  to  the  level 
of  approximation  that  is  involved  in  computing  vibrational  relaxation  and  its  effect  on 
dissociation. 
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II.  Mathematical  Formulation 

In  this  section,  we  introduce  the  set  of  coupled  partial  differential  equations  that  de¬ 
scribes  the  dynamics  of  the  flow-field.  For  the  high  enthalpy  flows  of  interest,  the  standard 
mass,  momentum  and  energy  equations  are  not  sufficient  to  describe  the  thennochemical 
state  of  the  gas  and  physics  involved.  The  high  temperatures  inherent  to  these  flows  trigger 
a  number  of  chemical  and  energy  relaxation  processes,  or  “real-gas”  effects,  and  the  flows 
arc  frequently  in  substantia!  nonequilibrium.  We  account  for  real- gas  effects  by  introducing 
finite- rate  chemical  reactions  and  internal  energy  relaxation  to  the  conservation  equations. 

The  Navier-Stokes  equation  set  has  been  extended  so  that  it  accurately  describes  the 
physics  of  a  gas  that  is  vibrationally  excited  and  chemically  reacting  [Candler  (1988), 
Olejniczak  (1997)],  and  it  is  only  briefly  presented  here.  The  equation  of  state  of  the  gas 
is  developed.  The  models  used  for  the  transfer  of  mass  between  different  chemical  species, 
transport  of  momentum,  and  the  transfer  of  energy  between  different  energy  inodes  are 
discussed  along  with  their  underlying  assumptions. 

1.  Basic  Assumptions 

The  flowfields  are  assumed  to  be  accurately  described  by  a  continuum  formulation. 
For  this  assumption  to  be  valid,  the  Knudsen  number  defined  as  the  ratio  of  the  mean  free 
path  to  the  length  scale  of  the  flow  Kn  —  X/L,  must  be  small.  For  the  flows  of  interest  the 
Knudsen  number  is  of  the  order  of  0.1-0.001  based  on  the  model  size,  and  thus  the  flowfield 
is  in  the  continuum  regime.  It  is  implicit  in  the  continuum  formulation  that  there  is  little 
statistical  variation  at  any  point  in  the  flow  and  therefore  the  continuum  description  of 
the  viscous  fluxes  is  consistent. 

We  restrict  the  formulation  to  non-ionized,  chemically  reacting  flows.  The  rotational 
state  of  the  gas  typically  equilibrates  within  a  few  collisions  and  we  can  assume  that  the 
translational  and  rotational  states  of  the  gas  are  in  equilibrium  at  a  common  temperature. 
In  this  work  we  do  not  consider  electronic  excitation  of  the  particles. 

2.  The  Conservation  Equations 

For  a  chemically  reacting  mixture  of  ns  species,  of  which  nd  are  diatomic,  there  are 
ns  species  mass  conservation  equations  of  the  form 

^  +  V  •  {psu)  =  -V  ■  ( psvs )  +  ws  , 

where  ps  is  the  density  of  species  s,  va  is  the  species  diffusion  velocity,  and  wa  is  the  rate  of 
production  of  the  species  due  to  chemical  reactions.  By  summing  the  species  density  over 
all  species  we  get  the  mixture  density  p.  If  we  sum  the  species  mass  conservation  equation 
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over  all  species  s,  the  source  terms  that  appear  on  the  right  hand  side  sum  identically 
to  zero  and  we  recover  the  traditional  continuity  equation.  The  models  used  to  describe 
the  diffusion  velocities  and  the  form  of  the  source  terms  will  be  discussed  in  subsequent 
sections. 

The  conservation  of  linear  momentum  is  expressed  as 

ot 

where  p  is  the  thermodynamic  pressure  and  r  is  the  viscous  stress  tensor.  This  vector 
equation  represents  3  momentum  equations,  one  for  each  spatial  dimensions  of  the  prob¬ 
lem.  The  above  expression  assumes  that  there  are  no  external  body  forces  exerted  on 
the  particles.  Note  that  diffusion  velocities  do  not  appear  in  the  momentum  equation  as 
their  contributions  sum  identically  to  zero  (see  Olejniczak  (1997)  for  the  details  of  the 
calculation). 

The  conservation  of  total  energy  of  the  mixture  is  given  by 

dE  ns^ 

—  +  V  •  {{E  +  p)u)  ~  V  ■  (fu)  -  V  ■  (gt  +  gr  +  qv)  -  V  •  J](psMs)  , 

S—l 

where  qt ,  qr,  qv  are  the  translational,  rotational  and  vibrational  heat  flux  vectors,  and  hs 
is  the  total  enthalpy  per  unit  mass  of  species  s.  The  terms  present  in  the  energy  equation 
will  be  discussed  separately. 

In  our  formulation  we  assume  that  vibration-vibration  energy  coupling  is  very  strong 
such  that  all  vibrational  modes  are  in  equilibrium  with  each  other.  Thus  the  vibrational 
energy  of  the  diatomic  species  can  be  described  by  a  single  vibrational  temperature.  This 
is  a  significant  simplification  as  we  only  have  one  vibrational  energy  conservation  equation, 
while  in  principle  we  need  to  have  a  separate  equation  for  each  of  the  polyatomic  species. 
The  equation  for  the  vibrational  energy  per  unit  volume  is  given  as 

0E 

+  V  •  (Evu)  =  -  V  ■  qv  -  V  •  ^(Pie^iTs)  +  wv  , 

S— 1 

where  eVs  is  the  vibrational  energy  per  unit  mass  of  species  s.  The  source  term  appearing 
on  the  right  hand  side  of  the  equation  will  be  discussed  in  the  subsequent  sections, 

3.  Equations  of  State 

In  this  section,  we  show  how  conserved  variables  relate  to  the  primitive  flow  variables, 
such  as  temperature,  and  pressure  that  appears  in  the  conservation  equations.  The  total 
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energy  per  unit  volume  of  the  gas  is  defined  as 


ns  j  nd  ns 

E  -  Y^PsCvaT  +  -p|u|2  +  Y^evsps  +  Ylh°Ps  ' 

5—1  5=1  5=1 


where  Cvs  is  the  translational-rotational  specific  heat  at  constant  volume  for  species  s, 
and  is  given  by 


cV3  =  aitr)+cv 


{rot) 

'5 


The  translational  and  rotational  specific  heats  are  given  by 


(tr)  _3  R 
Vs  “  2  Ms 

r<  (rot)  _ R 

°Vs  ~Ml 


s  <  nd  , 


where  R  is  the  universal  gas  constant  and  Ms  is  the  species  molecular  weight.  Thus  the 
first  term  in  the  energy  corresponds  to  the  internal  energy  of  the  gas  due  to  thermal 
translational-rotational  motion  of  particles.  The  second  term  represents  the  bulk  kinetic 
energy  of  the  gas.  The  third  term  represents  the  vibrational  energy  of  the  polyatomic 
species.  The  vibrational  energy  per  unit  mass  evs  of  species  s  is  related  to  the  vibrational 
temperature  Tv  assuming  a  Boltzmann  distribution  exists.  In  this  work,  we  use  a  simple 
harmonic  oscillator  model  to  represent  the  vibrational  excitation  of  diatomic  species  anti 
this  internal  energy  is  expressed  as 


R  dvs 

6va  ~  Ws  exp{8vJTv)-l  ' 

where  0V3  is  the  characteristic  temperature  of  vibration  for  species  s.  Note  that  it  is  not 
possible  to  solve  this  equation  directly  for  Tv,  and  therefore  the  vibrational  temperature 
must  be  computed  iteratively  for  a  given  vibrational  energy  Ev  =  ^s=1  psevs.  The  final 
term  in  the  expression  for  the  total  energy  represents  the  chemical  energy  of  formation 
associated  with  each  species. 


The  thermodynamic  pressure  p  is  related  to  the  translational- rotational  temperature 
through  a  perfect  gas  law  and  is  obtained  by  Dalton’s  Law  of  partial  pressures  of  the 
species  as 


rial  fib 

p=J2p*  =  Y.p>wT  =  pkT' 

S—  1  3=1  s 


and  R  =  ^  — 


ns  n 

_  ps  R 
P  Ms 


S  =  1 


The  enthalpy  per  unit  mass  of  species  s  is  defined  to  be 


4.  The  Diffusive  Terms 

In  this  section  we  describe  the  viscous,  thermal  and  mass  diffusion  fluxes  that  appear 
in  the  conservation  equations,  and  their  constitutive  relations. 

The  viscous  stress  tensor  is  formed  assuming  a  Newtonian  fluid  and  is  given  by 

f  =  p(Vu  +  (Vu)r)  +  A(V  ■  u)I  , 
with  the  Stokes  hypothesis  for  the  bulk  viscosity 

2 fi  +  3A  —  0  , 

where  p  is  the  kinematic  viscosity  of  the  mixture.  The  viscosity  / 1  depends  on  the  state  of 
the  fluid  and  its  composition  and  will  be  discussed  in  detail  below. 

The  heat  fluxes  are  assumed  to  obey  Fourier’s  Law  and  given  by 

9i  =  ~ '/CjVTi  , 

where  i  represents  any  of  the  internal  energy  modes  of  the  gas,  with  /t,  and  Ti  the  cor¬ 
responding  conductivity  and  temperature.  Since  the  translational  and  rotational  states 
of  the  gas  are  in  equilibrium  at  a  common  temperature,  it  is  convenient  to  combine  the 
corresponding  heat  fluxes  into  a  sum  of  the  conductivities  times  the  temperature  gradient 
as 

q  =  qt+qr  =  -(/c(  4-  Kr)  VT  . 

The  thermal  conductivities  for  the  translational,  rotational,  and  vibrational  energy  modes 
are  determined  from  an  Eucken  relation  [Vincenti  and  Kruger  (1986)]  as 

-  §A**C’v?r)  .  =  HsC^ot)  ,  «vS  =  »sCv<?ib)  , 

where  fis  and  ks  are  the  species  viscosities  and  thermal  conductivities.  The  vibrational 
specific  heat  is  computed  directly  using 

/-i  (uift)  _ 

vs  ~  dTv  ' 

Mass  diffusion  is  driven  by  gradients  of  concentration,  pressure,  and  temperature. 
However,  for  the  flows  of  interest  only  the  terms  due  to  concentration  gradients  are  sig¬ 
nificant,  Thermal  diffusion  is  important  in  regions  where  gradients  are  large,  but  in  these 
regions  the  fundamental  continuum  assumption  is  typically  invalid  and  we  can  no  longer 
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apply  the  Navier-Stokes  equations  [Boyd  et  al.  ( 1995) j .  The  baro-diffusion  terms  are  negli¬ 
gible.  Therefore,  for  this  work  the  mass  diffusion  is  solely  based  on  concentration  gradients 
and  Pick’s  Law  is  used.  The  mass  diffusion  flux  is 


psva  =  -pVaV{—)  . 

P 


The  multicomponent  diffusion  coefficients  Vs  are  replaced  with  a  single  binary  diffusion 
coefficient  V,  which  is  derived  assuming  a  constant  Lewis  number,  Le,  defined  as 


Le  =  D 


pCv 


K 


where  Cp  is  the  translational- rotational  specific  heat  at  constant  pressure  and  k  is  the 
thermal  conductivity  of  the  mixture.  The  effects  of  multicomponent  diffusion  have  been 
neglected,  however  they  should  be  included  if  light  species,  such  as  hydrogen,  are  present 
in  the  flow. 


It  now  remains  to  provide  expressions  for  the  mixture  viscosity,  the  thermal  conduc¬ 
tivity  and  the  species  viscosities.  The  mixture  viscosity  p  and  thermal  conductivity  k  are 
obtained  from  the  species  viscosities  and  conductivities  fia>  ns  using  Wilke’s  semi-empirical 
mixing  rule  [Wilke  (1950)] 


*-E 


Xsf-ts 


«-£ 


<Ps 


where 


r 


j 


and  cs  =  psj  p  is  the  mass  fraction  of  species  s.  The  model  of  Biottner  is  used  for  the 
species  viscosities  p3,  which  are  calculated  from  curve  fits  of  the  form 


ps  =  0.1  exp[(As  In  r  +  Bs)  In  T  +  Cs]  , 

where  As,  Bs,  and  Ca  are  constants.  This  model  is  valid  for  a  large  range  of  temperature 
(to  about  10,000  K),  and  is  suitable  for  this  work. 


5.  The  Chemical  Source  Terms 

Each  of  the  source  terms  ws  that  appear  in  the  mass  conservation  equations  represents 
the  rate  of  production  and  destruction  of  species  s  due  to  chemical  reactions.  In  this  section 
we  present  the  chemical  source  terms  that  appear  in  the  mass  conservation  equations. 
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In  the  high  enthalpy  flows  that  we  examine  in  this  work,  two  different  gases  are  used: 
nitrogen  and  air.  The  treatment  of  nitrogen  can  be  considered  as  a  special  case  of  the 
air  gas,  and  thus  the  set  of  conservation  equations  for  nitrogen  flows  is  a  subset  of  the 
equation  set  that  describes  air  flows.  For  high  temperature  air  with  no  ionization,  there 
are  five  chemical  reactions  that  are  important,  involving  five  species:  N2,  02,  NO,  N,  and 
0.  The  chemical  reactions  involved  are 


N2  +  M  2N  +  M 
02  +  M  20  +  M 


NO  +  Mf^N  +  0  +  M 


N2+O^NO  +  N 
NO  +  O  02  +  N 


The  reactions  are  presented  so  that  they  are  endothermic  in  the  forward  direction.  The 
first  three  reactions  represent  dissociation  of  the  diatomic  species  N2,  02,  and  NO  due 
to  collisions.  In  these  reactions,  M  is  arbitrary  and  represents  the  collision  partner  which 
provides  the  energy  to  break  the  chemical  bond.  The  three  corresponding  reverse  reactions 
are  three-body  recombination,  during  which  the  collision  partner  M  absorbs  the  energy 
released  during  the  formation  of  chemical  bonds.  The  two  remaining  reactions  are  the 
Zeldovich  exchange  reactions,  which  are  the  primary  source  of  nitric  oxide  (NO)  in  a 
hypersonic  flow. 

The  rates  of  the  forward  and  backward  reactions  take  an  Arrhenius  form  governed  by 
the  forward  and  backward  rate  coefficients  kjVm  and  kbrTn  respectively.  We  may  write  the 
rate  of  each  reaction  as  a  sum  of  the  forward  and  backward  rates 
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The  chemical  source  terms  can  now  be  expressed  in  terms  of  the  individual  reaction  rates 

=  MnjO^i  +  ^4) 

™o2  =  K>20^2  ~Tlb) 
wNO  =  Mho(Hz  -T^4  +^5) 
irjN  =  MN(—2Tl\  —  7ls  —  IZ4  —  7^5) 
yj0  =  M0  (—27^2  —  7^3  -j-  TL\  +  7^) 

We  note  that  the  sum  of  all  the  source  terms  is  identically  zero,  so  that  conservation  of 
total  mass  in  satisfied,  and  elemental  conservation  of  atomic  nitrogen  and  oxygen  is  also 
satisfied. 

The  forward  and  backward  reaction  rate  coefficients  are  in  general  affected  by  the 
level  of  thermal  nonequilibrium  in  the  flow.  In  this  work  we  use  the  model  of  Park  (1986) 
for  vibration-dissociation  coupling.  The  Park  model  has  been  used  extensively,  since  it 
is  simple  to  implement  and  has  been  shown  to  give  satisfactory  results  for  a  variety  of 
How  conditions  [Olejniczak  (1997)].  The  Park  model  assumes  that  the  reaction  rates  are 
functions  of  an  effective  temperature,  and  therefore  we  represent  the  forward  reaction  rate 
coefficients  by  the  modified  Arrhenius  form 

kSm  =  CfmVfi  eXP(-0m/Ten)  , 

where  C/m,  rfm,  and  dm  are  found  from  curve  fits  to  experimental  data.  In  this  work,  the 
values  recommended  by  Park  (1988)  are  used.  The  effective  temperature  T ^  that  appears 
in  the  expression  is  some  function  of  the  translational-rotational  and  vibrational  temper¬ 
atures.  Park  recommends  that  the  geometric  average  of  the  translational  and  vibrational 
temperatures  is  used  for  the  forward  rates  of  the  dissociation  reactions,  as 

Teff  =  \/TTv  . 


The  backward  reaction  rate  coefficients  are  related  to  the  forward  reaction  rate  coefficients 
through  the  equilibrium  constant  defined  as  their  ratio 


K*JT)  = 


The  equilibrium  constant,  Kcq,  is  also  a  curve  fit  to  experimental  data,  and  is  a  function 
only  of  T 

=  Cmexp(Aim  +  A^mZ  +  AamZ2  +  A4mZ3  +  AsmZ4)  , 

where  Z  =  10,000 /T.  The  constants  that  appear  in  the  expression  for  K<, q  are  given  by 
Park  (1985). 
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6.  Vibrational  Energy  Source  Term 

The  source  term  that  appears  on  the  right  hand  side  of  the  vibrational  energy  equation 
(a;,;)  is  presented  here.  This  source  term  is  due  to  the  exchange  of  energy  between  the 
translational- rotational  and  vibrational  modes  of  the  gas.  Additionally,  the  creation  of 
diatomic  molecules,  due  to  recombination  reactions,  contributes  to  the  vibrational  energy 
source  term. 

There  are  many  energy  exchange  mechanisms  between  the  different  internal  energy 
modes  of  the  gas,  including  exchange  between  translational  and  vibrational  modes,  ex¬ 
change  between  rotational  and  vibrational  modes,  and  exchange  between  vibrational  modes 
of  the  separate  species.  In  this  work,  we  have  assumed  that  the  vibration-vibration  ex¬ 
changes  are  very  fast,  and  thus  the  vibrational  modes  of  the  gas  are  in  equilibrium  at  the 
same  vibrational  temperature  TV  The  vibration-translation  and  vibration- rotation  energy 
exchange  rates  are  combined  in  a  single  energy  exchange  rate  Qt~v,  and  the  Landau- Teller 
model  is  used  [Vincenti  and  Kruger  (1986)].  The  energy  exchange  rate  is  given  as 


ALT) 


Qt-v  =Ps- 


6t/i 


<TS> 


where  e’s  is  the  vibrational  energy  per  unit  mass  of  species  s  evaluated  at  the  translational 
temperature  (equilibrium),  and  <  t3  >  is  the  molar  averaged  Landau-Teller  relaxation 
time,  and  is  given  by  Lee  (1985)  as 


<  Ts  >  = 


T.  rXT 

£.  Xr/r„  ’ 


where  Xr  has  been  defined  previously  and  Tsr  is  Landau-Teller  inter-species  relaxation 
time  given  by  Millikan  and  White  (1963),  as 

Tsr  ~  -  exp[Asr(T~1/f3  -  0.01  5/j^4)  -  18.42]  ,  (p  in  atm) 

P 

Aar  =  1.16  x  10 -3pl'20vi/3, 


Psr  = 


M3  Mr 
Ms  +  Mr 


Part  of  the  vibrational  energy  source  terms  is  due  to  the  creation  and  destruction  of 
diatomic  species  in  chemical  reactions.  For  this  work,  we  assume  that  the  molecules  that 
are  recombining  and  molecules  that  dissociate  carry  vibrational  energy  given  by  the  local 
vibrational  temperature.  Therefore,  the  rate  of  production  of  vibrational  energy  due  to 
chemical  reactions  is  just 

Qva  =  V}aevs  . 
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The  vibrational  energy  source  term  ws  is  simply  the  sum  of  the  contributions  of  the  two 
processes,  namely  recombination  and  translation- vibration  energy  exchange,  and  is  given 
by 

nd  nd 

E  <c,T) 

Qt-v  $  ■ 

3—1  5=1 
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III.  Numerical  Method 


In  this  section  we  present  the  numerical  method  we  use  to  solve  the  coupled  partial 
differential  equations  that  describe  the  flows  of  interest.  The  conservation  laws  presented 
in  the  previous  section  were  expressed  in  their  strong  form.  However,  we  obtain  numerical 
solutions  to  this  set  of  conservation  equations  by  considering  the  weak  form  of  the  equa¬ 
tions.  VVe  use  a  standard  finite  volume  method  to  discretize  the  governing  equations  in 
the  approximation  space. 

We  introduce  the  finite  volume  formulation,  and  the  methods  we  use  to  evaluate  the 
inviscid  and  viscous  terms.  The  implicit  method  employed  for  time  advancement  is  also 
briefly  presented.  We  will  restrict  the  remaining  of  this  section  to  the  two-dimensional 
case  in  order  to  simplify  the  calculations,  but  it  should  be  noted  that  it  is  straightforward 
to  extend  the  numerical  method  to  three  dimensions. 


I  .  The  Finite  Volume  Method 

In  this  section,  we  introduce  the  finite  volume  method  of  discretization.  We  begin 
with  the  strong  form  of  the  equations  presented  in  the  previous  section.  We  write  the 
vector  of  conserved  variables  as 

U  =  (/>!,  -  -  -  ,P™,pu,pv,Ev,E)T  , 


and  the  conservation  equations  can  now  take  the  compact  form 

dU 


dt 


+  V-F  =  W  , 


where  F  is  the  corresponding  total  Hux  vector  and  W  is  the  vector  of  source  terms 

W  =  (u>i, . . . ,  rcns,0,0,  u?d,  0)T  . 

The  total  flux  can  be  split  into  the  convective  (inviscid)  and  diffusive  flux  vectors 

F  =  Fi  +  Fv  . 


The  Cartesian  components  of  the  convective  and  diffusive  flux  vectors  are 


/  p\u  \ 

Piv  \ 

p2U 

P2V 

,  Gi  = 

PnsV 

pun  +  p 

puv 

pvu 

pvv  +  p 

EylL 

Evv 

\  ( E  +  p)uj 

\  (E  +p)vj 

18 


F„  = 


P  lUl 
P2U2 

AisUns 

^TT 


—  T 


Ey 


9vx  +  Es=x(^et>sus) 

V  Qx  +  9ujc  —  (Tti)z  +  (Ps^sU-j)  / 


/ 


Gv  — 


P 1V1 
P2V2 


Pnsvns 


-T, 


Qvy  +  Z^S=l(p5ev^Vs) 

\9v  +  Qvy  -  (fu)y  +  i(/>jA*v,)  / 


yy 


In  the  finite  volume  formulation  the  weak  form  of  the  conservation  equations  is  ob¬ 
tained  by  integrating  over  an  arbitrary  control  volume  Q,  which,  for  the  purposes  of  this 
work,  is  assumed  to  be  fixed  in  space.  Green's  Theorem  is  applied  to  convert  the  volume 
integral  containing  the  flux  vector,  to  a  surface  integral  over  the  surface  that  encloses  SI 
The  resulting  weak  form  of  the  equation  is 

w  +  LfjP.W-*, 

where  V  is  the  total  volume,  n  is  the  outward-pointing  normal  to  the  surface,  and  U,  W 
are  averaged  over  fl.  We  discretize  physical  space  by  dividing  it  into  a  series  of  volume 
elements,  which  can  be  arbitrary  polygons.  We  do  this  by  building  a  regular  mesh  over 
the  domain  of  interest,  resulting  in  a  structured  grid.  The  grid  directions  are  defined  such 
that  i  index  increases  in  the  ^-direction  and  j  increases  in  the  positive  redirection.  Each 
of  the  volume  elements  in  the  two-dimensional  domain  is  a  quadrilateral,  and  represents  a 
cell  in  computational  space.  Flow  variables  are  stored  at  the  cell  centers  rather  than  the 
nodes.  Quantities  at  the  cell  faces  are  computed  from  the  solution  reconstruction  based 
on  the  cell-centered  data.  The  time  rate  of  change  of  0  in  the  2-D  volume  element  ij  is 
represented  as  a  sum  of  the  fluxes  over  the  four  cell  faces  and  the  generation  of  U  within 
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the  volume 


BOij  \ 

dt 


=  - y-  [(?  ■  ft)i+i  ,j  -  (F  ■ 

+  {F-h)iJ+iSitj+i  -(F  -n)itj_ iSij-} 


+  Wi 


j.j 


The  ±1/2  indices  on  the  flux  and  surface  area  terms  indicate  that  quantities  should  be 
evaluated  at  the  appropriate  cell  face.  This  expression  is  the  discretized  form  of  the 
conservation  equations  for  a  two-dimensional  domain.  The  extension  of  the  formulation  to 
three  dimensions,  where  the  volume  elements  are  hexahcdra,  is  straight-forward  and  results 
in  an  analogous  expression.  The  axisymmetric  formulation,  however,  is  non- trivial,  and 
ultimately  yields  an  expression  similar  to  that  of  the  two-dimensional  case.  Integrating  the 
conservation  equations  over  an  axisymmetric  domain  introduces  extra  terms  that  appear 
as  sources  in  the  forcing  vector  W . 


2.  Evaluation  of  the  Fluxes 

In  this  section  we  present  the  methods  we  use  for  evaluating  the  fluxes  that  appear 
on  the  right  hand  side  of  the  discretized  form  of  the  conservation  equations.  We  have 
decomposed  the  fluxes  into  the  convective  part  (which  consists  of  the  advection  terms 
and  the  pressure)  and  diffusive  part.  This  decomposition  is  useful  because  we  treat  the 
convective  and  diffusive  fluxes  differently.  The  in  vise  id  terms  make  the  system  of  equations 
hyperbolic,  while  the  viscous  terms  are  diffusive  in  nature,  and  make  the  system  parabolic 
when  combined  with  the  transient  term.  This  complicates  the  evaluation  of  the  right  hand 
side  of  the  conservation  equations.  The  elliptic  nature  of  the  diffusive  terms  allows  for  a 
simple  central  difference  scheme  to  be  used  [Hirsch  (1991)].  A  more  sophisticated  approacli 
is  necessary  for  the  calculation  of  the  inviscid  terms,  in  order  to  maintain  stability  of  the 
numerical  scheme. 

In  the  computation  of  the  inviscid  terms,  we  make  use  of  the  fact  that  in  hyperbolic 
systems  of  equations,  information  propagates  along  definable  characteristic  directions  in 
space-time.  This  is  commonly  referred  to  as  upwind  biasing.  We  employ  a  modified 
form  of  the  Steger- Warming  (1981)  flux-vector  splitting,  a  method  widely  used  due  to  its 
robustness,  which  is  based  on  characteristic  theory.  We  present  the  application  of  this 
method  specifically  to  the  set  of  conservation  equations  relevant  to  this  work  (see  also 
Wright  (1997)),  but  it  should  be  noted  that  the  method  has  been  implemented  for  more 
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complex  sets  of  conservation  equations.  We  consider  the  in  viscid  dux  through  a  face 


Pi«'  \ 

P2U* 

pnsu' 

puu'  +  ps'x 
pvu'  +  ps'y 
Evu' 

(E  +  p)v!  / 

wliere  are  the  components  of  the  unit  vector  normal  to  the  surface  in  question,  and 
v!  is  the  component  of  the  velocity  normal  to  the  surface.  We  observe  that  the  flux  F\  is 
linear  and  homogeneous  in  the  vector  of  conserved  variables  U,  i.e., 

FKW)  =  XFj(U)  , 

where  X  is  an  arbitrary  scalar.  This  can  be  shown  if  we  express  the  iti viscid  flux  F'{  in 
terms  of  the  components  of  U .  The  homogeneity  of  the  flux  allows  as  to  express  it  exactly 
in  terms  of  the  linearization 

f)  Ff 

F’,(U)  =  -gjjU  =  flU  , 

where  A!  is  the  inviscid  flux  Jacobian  matrix  in  the  direction  of  h.  We  follow  the  approach 

of  Steger  and  Warming  to  split  the  inviscid  fluxes  into  positively- moving  and  negatively- 

moving  components  using  an  upwind  biasing  scheme.  The  partitioning  is  performed  by 

diagonalizing  A',  and  the  fluxes  are  split  according  to  the  sign  of  the  eigenvalues  of  the 

dF'f 

Jacobian.  The  diagonalization  of  is  rather  cumbersome,  but  we  can  simplify  the 
operation  by  breaking  the  Jacobian  into  components,  as 

,  _  dU  dV  dFj  dV 
~  dV  dU  dV  d~U  ’ 

where  V  is  a  vector  of  primitive  variables,  introduced  for  convenience.  The  choice  of  V  is 
not  unique,  but  for  the  particular  inviscid  fluxes  it  is  convenient  to  use 

^  (pi  i  Pi i  ■  ■  ■  j  Pnst  ev ,  p)  , 


Fi  =  Ft  n  = 


where  ev  is  the  total  vibrational  energy  per  unit  mass.  The  transformation  matrices 
between  primitive  and  conservative  variables  are  denoted  by  5_1  and  S.  We  should,  in 
principle,  be  able  to  diagonalize  using  a  similarity  transformation 


dV  dFi 
dU  dV 


CJAa,Ca>  , 
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where  is  a  diagonal  matrix  of  eigenvalues,  C^}  is  the  matrix  of  eigenvectors  of  A\ 
and  Ca'  such  that  C^}CA‘  —  /•  The  diagonal  izat  ion  process  has  now  been  simplified, 
however  the  resulting  expressions  for  the  eigenvectors,  although  mathematically  correct, 
may  exhibit  singularities  and  are  not  suitable  for  use  in  a  numerical  algorithm.  Therefore, 
we  introduce  an  orthogonal  transformation  R,  which  transforms  the  momenta  from  a 
coordinate  system  locally  aligned  with  the  face  in  question,  to  the  Cartesian  system.  We 
can  now  write  the  Jacobian  matrix  as 

A1  =  S~l R-'C^1  Aa'CaR  S  , 


where  Ca  does  not  exhibit  singularities.  The  diagonal  components  of  AA'  are  the  con¬ 
vection  speeds  of  the  characteristic  variables.  The  flux  vector  can  now  be  split  into  the 
positively  and  negatively  moving  components  defined  by 

f;  =  s~{  r-'c^1  a+car  su  =  a'+u 

FL  =  S~x R~lC^x A.CaR  S  U  -  A'_U  , 

where  A+  and  A  are  the  split  eigenvalue  matrices  consisting  of  the  positive  and  negative 
eigenvalues  respectively.  The  total  inviscid  Hux  vector  through  the  face  in  question  can  be 
expressed  as  the  sum  of  positively  and  negatively  moving  components 


f;  =  F|  +  FL  . 


In  this  basic  discretization,  the  split  fluxes  at  each  cell  face  are  evaluated  using 


ni+w 


=  A+tjViJ  * 
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where  A'±  and  U  are  evaluated  based  on  the  upwind  cell  data,  MacCormack  and  Candler 
(1989)  proposed  a  modification  to  this  method,  which  reduces  the  amount  of  numerical 
dissipation  of  the  original  method.  Their  modification  suggests  that  the  Jacobian  matrix 
should  be  evaluated  at  the  cell  face,  rather  than  using  the  upwind  cell  data.  This  approach 
works  well  in  regions  of  weak  gradients,  but  additional  dissipation  is  required  to  capture 
strong  gradients  and  shock  waves.  In  this  work,  we  use  a  hybrid  method,  in  which  we 
use  pressure  weighted  average  quantities  between  adjacent  cells  to  evaluate  the  Jacobian. 
The  method  smoothly  switches  from  modified  to  true  Steger- Warming  in  regions  of  large 
pressure  gradients.  This  flux  evaluation  method  as  originally  formulated  is  only  first  order 
accurate  in  space,  however  it  is  possible  to  achieve  higher  order  of  accuracy.  Higher  order 
accurate  schemes  can  produce  more  resolved  solutions  on  a  given  computational  mesh.  We 
use  up  to  second  order  accurate  fluxes  in  this  work.  We  obtain  higher  order  of  accuracy  by 
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extrapolating  the  conserved  variables  to  the  cell  face  based  on  the  purely  upwind  cell  data. 
We  detect  the  presence  of  shock  waves  by  identifying  pressure  jumps  greater  than  50%  of 
the  smallest  pressure  involved  in  the  flux  evaluation  stencil,  in  which  case  we  switch  the 
stencil  to  first  order  For  the  details  of  this  formulation  see  Druguct  et  ai  (2003). 

If  we  substitute  the  split  fluxes  into  the  discretized  form  of  the  conservation  equations, 
we  get  the  upwind  finite  volume  formulation  of  the  in  viscid  equations 


dUi,j 

dt 


-W-i-ijSi-uUu  -  A'-i+iJS«i 

+  +  $Ui,3  ~  B+i  j-%  $ i,j  ~  3  ' “  !  ) 

«U  -  BLiJHSiJ+iUw)]  +  , 


where  we  have  implicitly  defined  A*  and  B*  to  denote  the  in  viscid  flux  Jacobi  ans  in  the  £ 
and  tj  directions  respectively. 

The  elliptic  nature  of  the  viscous  fluxes  makes  their  evaluation  possible  with  a  simple 
finite  difference  scheme.  The  viscous  flux  vector  is  evaluated  at  the  required  cell  faces 
using  central  differencing.  The  viscous  fluxes  integrated  over  a  face  take  the  simple  form 

"  (i^  *  j^^Si j+i  , 

The  contribution  of  the  viscous  fluxes  is  summed  over  the  four  cell  faces  and  added  to  the 
right  hand  side  of  the  discretized  form  of  the  conservation  equations.  The  required  deriva¬ 
tives  that  appear  in  the  viscous  fluxes  are  taken  with  respect  to  the  general  curvilinear 
coordinate  system  and  are  then  transformed  to  the  Cartesian  system.  The  calculation 
of  these  terms  can  be  found  in  Hirsch  (1991), 


3*  Time  Advancement 

In  this  section  we  discuss  how  the  solution  is  integrated  in  time  to  its  steady  state. 
We  are  interested  solely  in  the  steady  state  solution,  if  one  exists,  and  thus  time  accuracy 
is  not  important  for  our  computations.  Therefore,  we  employ  an  implicit  method,  which 
allows  for  large  time  steps  to  be  taken,  while  maintaining  stability  of  the  numerical  scheme. 
We  give  a  short  description  of  the  implicit  formulation  here.  We  first  introduce  the  explicit 
method  of  solution,  and  then  we  extend  it  to  the  implicit  formulation. 

We  discretize  the  transient  term  that  appears  in  the  conservation  equations,  in  order  to 
time-march  the  solution  toward  the  steady  state.  We  use  simple  first-order  Euler  backward 
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differencing,  in  which  the  time  derivative  becomes 


dUij 

at 


u. 


n  +  1 


A  Un* 

iJ-  +  0{At2)  ~ 


At 


and  the  superscript  denotes  the  time  level  in  the  temporal  discretization.  This  expression 
is  substituted  directly  into  the  discretized  form  of  the  conservation  equations,  and  the 
solution  is  marched  forward  in  time  by  discrete  time  steps.  The  explicit  formulation  of  the 
inviscid  problem  is 


AU*,  = 


-Vz[(A'+<H.iS<+i'iu‘J  -  +i-i 

~  (A'~  i -  \ ,jSi~  3 ,3Ui>J  ~  5i+  3-J  ^i+ 1  ) 


where  AU?j  is  the  explicit  change  in  the  solution  vector  (residual)  at  cell  ij  at  time  step 
n.  The  solution  is  updated  to  the  next  time  level  using 


Utf^U^  +  AUrj 


We  call  this  an  explicit  method  because  the  solution  in  the  next  time  level  n+1,  is  computed 
explicitly  by  evaluating  the  right  hand  side  at  the  current  time  level  n.  The  solution  at 
any  point  at  time  level  n  +  1  is  independent  of  the  solution  at  other  points  at  the  same 
time  level. 

With  a  fully  implicit  method,  the  time-step  limitations  inherent  to  explicit  solvers 
can  be  overcome.  The  use  of  an  implicit  solver  can  allow,  in  general,  for  time  steps  much 
larger  than  the  stable  explicit  time  step.  In  an  implicit  method,  the  solution  at  any  point 
depends  on  the  solution  of  all  other  points  at  the  time  level  n+1.  Therefore,  in  order  to 
employ  implicit  time  advancement,  it  is  necessary  to  evaluate  the  right  hand  side  of  the 
conservation  equations  at  time  level  n  +  1 .  We  express  the  fluxes  that  appear  on  the  right 
hand  side  in  terms  of  the  linearization 


f;”+i 


5f; 


dU 


(u 


n  +  1 


+  A'n6Un 


-Un)+0{At2) 


where  we  have  defined  the  vector  SUn  =  Un+]  -  Un  and  we  have  frozen  the  variation  of 
the  inviscid  flux  Jacobian  A,n  between  time  steps.  We  apply  this  linearization  to  the  split 


24 


fluxes  and  source  vector  to  obtain  the  fully  implicit  upwind  formulation  of  the  in  viscid 
problem 

-{A'-i-ijSi-ijWiJ  -  A'-<+ijs<H ,/«+>.;) 

-Ai  CTjitTj  =  A£/£  , 

where  C?j  is  the  Jacobian  of  W  with  respect  to  U.  In  this  implicit  equation,  A£/“j-  is  the 
explicit  residual  as  defined  previously. 

We  include  the  implicit  approximation  of  the  viscous  terms  in  the  implicit  equation 
by  applying  the  linearization  to  the  viscous  fluxes.  We  define  SF^n  and  8 G'vn  such  that 


F'n+1  =  F'"  +  SFV 

G'vn+l  =  c;n  +  5<?;n 


/n 


Therefore,  we  only  need  to  find  expressions  for  SF 'v  and  8G'V.  We  note  that  F'v  and  G'v 
are  functions  of  U  and  its  spatial  derivatives.  The  derivatives  are  taken  with  respect  to 
the  general  curvilinear  coordinate  system,  and  then  transformed  to  the  Cartesian  system 
using  the  operators 

d  _  d^  d  dr]  d  d  _  d£  d  dr]  d 

dx  dx  d£  ^  dx  drj  ’  dy  dy  d£  dy  dr) 

We  apply  the  thin-layer  assumption  to  the  derivatives  in  the  implicit  viscous  terms.  In 
this  way,  we  selectively  retain  derivatives  depending  on  which  viscous  flux  we  wish  to 
approximate.  For  example,  in  the  approximation  of  SG’V  we  compute  the  derivatives  by 

d  ^  dr]  d  d  ~  drj  d 

dx  ~  dx  dt]  ’  dy  ~  dy  drj 

Similarly,  derivatives  with  respect  to  £  are  retained  when  approximating  <5 F‘v .  If  wc  assume 
that  the  terms  not  involving  derivatives  are  locally  constant,  we  can  write  F'  and  G'v  as  a 
function  of  the  derivatives  of  the  flow  variables  with  respect  to  £  or  i]  only,  as 

SK  =  M(  1  (iV)  ,  SG'V  -  M,  A(iV)  , 

and  we  introduce  for  convenience  the  vector  of  non-conscrved  variables 


V  =  (ci,c2, .. .  ,cns,  u,v,ev,T)r  . 
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The  matrices  and  M,t  can  be  found  in  Candler  (1988).  We  ultimately  express  the 
implicit  viscous  fluxes  in  terms  of  6U  by  mapping  from  <5V  to  SU  using  the  Jacobian 
N  —  and  we  obtain  the  final  form 

6F'  =  M^(N5U)  ,  SG'V  =  M^NSU)  , 

We  should  emphasize  here  that  ail  of  the  viscous  derivatives  are  retained  in  the  computation 
of  F'vn  and  G'vn,  and  thus  all  of  the  viscous  terms  influence  the  converged  solution.  With 
this  treatment  of  the  viscous  terms,  the  implicit  equation  for  the  inviscid  problem  away 
from  solid  boundaries  will  be  unchanged  if  we  simply  replace  the  inviscid  Jacobians  A '  and 
B'  with  A  and  B,  where. 

A+  =  A+  -  M^N  =  A _  +  M^N 

B+  =  B+-  M„N  B.=B.  +  M,tN 

The  implicit  formulation  involves  the  solution  of  a  linear  system  of  equations,  which 
has  a  strong  diagonal  dominance  due  to  the  nature  of  the  Steger- Warming  Jacobians.  In 
principle,  we  can  solve  exactly  for  the  6U?j  by  a  direct  inversion  of  the  operator. 

4.  The  Data-Paralle!  Line- Relaxation 

In  this  section  we  briefly  present  the  method  of  solution  of  the  implicit  equation.  We 
use  the  Data-Parallel  Line  Relaxation  method  of  Wright  et  al.  (1998).  This  method  was 
developed  to  take  advantage  of  parallel  computing  machines,  and  appears  to  be  superior 
to  other  implicit  methods  [Wright  (1997)].  Only  a  brief  description  of  this  method  is  given 
here.  Details  on  the  performance  of  the  method  can  be  found  in  Wright  (1997)  and  Wright 
ef  al.  (1998). 

The  standard  implicit  upwind  formulation  of  the  Navier-Stokes  equations  for  a  reacting 
gas  is  essentially  a  linear  system  of  equations.  The  solution  of  such  system  can,  in  principle, 
be  obtained  directly  by  inversion  of  a  large  block  banded  matrix.  This  becomes  very 
computationally  intensive,  as  such  an  inversion  must  be  performed  at  every  time  step  until 
convergence  is  achieved.  Therefore,  most  implicit  methods  seek  to  make  simplifications  in 
order  to  make  the  solution  of  the  problem  less  expensive.  Such  simplifications  are  possible 
only  on  the  left-hand  side  of  the  equations,  which  resembles  the  numerics.  The  right-hand 
side  of  the  equation  (residual)  should  be  unaltered. 

It  is  convenient  to  introduce  some  additional  notation  in  order  to  group  terms  on  the 
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left-hand  side  of  the  implicit  equation.  We  let 


AlJ 

~  A t1  +  (A  +i+±>iSi+i'j 

-,4'\ 

1 

-B,n. 
~  1 

* J—  3  _  5 

We  can  simplify  the  operator  if  we  assume  that  the  physical  problem  is  strongly 
coupled  in  the  t]  direction.  Then,  we  can  move  some  of  the  off-diagonal  terms  to  the 
right-hand  side  of  the  implicit  equation,  keeping  only  constant  i  terms  on  the  left-hand 
side,  In  this  way,  it  is  possible  to  solve  for  all  j  points  at  each  i  location  as  a  series  of  fully 
coupled  block  tri-diagonal  systems  aligned  in  the  r ]  direction.  We  relax  the  off-diagonal 
terms  by  updating  the  unknowns  on  the  right-hand  side  of  the  equation  with  the  SUij  we 
have  just  computed.  We  do  not  need  to  obtain  an  exact  solution  to  the  implicit  equation, 
but  rather  to  maintain  stability,  and  thus  only  a  few  iterations  are  required  per  time  step. 
The  solution  procedure  is  as  follows.  First,  the  implicit  terms  on  the  right-hand  side 
are  neglected  and  the  resulting  block  tri-diagonal  system  is  factored  and  solved  for  6U ^ 
according  to: 


1  +  -  CuW.7-1  = 


r(0) 


(0)  _ 


where  the  superscript  denotes  the  sub-iteration  number  (zero  denotes  the  initialization). 
Then,  a  series  of  kmax  relaxation  steps  is  performed,  where  we  solve 


1  +  Mjtufl  -  Cid6Utj’_,  =  -DtjSU&y  +  BtjSUil Z1  +  Al/Tj 


TW 


tW  _ 


(k-l) 


r(k-\) 


and  the  superscript  k  -  1  indicates  that  data  from  the  previous  relaxation  step  is  used. 
Finally,  the  solution  is  obtained  by 

8UTtx  =  sulk™az)  ■ 


There  are  two  advantages  to  this  method  of  solution.  First,  the  method  can  be  imple¬ 
mented  efficiently  on  a  distributed  memory  parallel  computer.  This  is  because  we  can  split 
the  problem  among  processors  in  the  i  direction,  effectively  assigning  a  series  of  i-columns 
to  each  processor.  Therefore,  each  relaxation  step  can  be  performed  simultaneously  as 
there  are  no  data  dependencies,  while  the  boundary  data  can  be  efficiently  transfered. 
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Secondly,  there  is  no  bias  in  the  solution  procedure,  which  has  been  shown  to  impose  a 
bias  in  the  solution  itself. 

5.  Boundary  Conditions 

In  this  section  we  describe  the  boundary  conditions  used  for  both  the  explicit  residual 
and  the  implicit  operator,  and  how  they  are  implemented  in  the  numerical  algorithm. 

In  the  finite  volume  method,  explicit  boundary  conditions  are  imposed  using  “dummy” 
cells.  An  extra  layer  of  cells  is  included  around  the  grid  that  represents  physical  space,  in 
which  variables  are  stored  in  order  to  impose  the  appropriate  boundary  conditions  for  the 
problem.  The  dummy  cells  are  constructed  in  such  a  way  that  derivatives  taken  in  terms 
of  the  curvilinear  body-fitted  coordinate  system  are  consistent  with  the  type  of  boundary 
condition  of  the  viscous  fluxes.  Because  we  split  the  fluxes  into  the  viscous  and  inviscid 
parts,  we  treat  the  boundary  conditions  separately. 

In  the  explicit  formulation,  we  simulate  no-slip  conditions  by  setting  the  velocity  at 
the  wall  to  zero,  and  extrapolating  the  velocity  components  at  the  boundary  cell  i,  1  based 
on  the  values  at  the  cell  i, 2.  We  simulate  isothermal  wall  conditions  in  a  similar  fashion. 
The  slip  conditions  require  that  we  have  a  model  for  how  the  velocity  component  tangent 
to  the  wait  is  computed  based  on  the  interior  values.  The  details  of  the  model  employed 
were  presented  in  the  previous  section.  Once  the  slip  velocity  has  been  determined,  we 
extrapolate  values  to  the  dummy  cells,  as  we  did  for  the  no-slip  case.  We  impose  internal 
energy  jumps  at  the  wall  in  exactly  the  same  fashion.  The  non-catalycity  of  the  wall  is 
imposed  by  setting  the  values  of  the  species  densities  at  the  dummy  cells  equal  to  the 
interior  values,  effectively  imposing  a  zero  mass- fraction  gradient.  Similarly,  fully  catalytic 
boundary  conditions  are  imposed  be  extrapolating  mass-fractions  to  the  dummy  cell  as¬ 
suming  the  mass- fractions  of  atomic  species  is  zero  at  the  wall.  Subsequently,  the  mass 
flux  of  polyatomic  species  at  the  wall  is  corrected  such  that  we  have  elemental  conservation 
during  this  process. 

The  implicit  treatment  of  the  boundary  conditions  is  analogous  to  that  of  the  explicit 
method.  In  the  implicit  formulation,  however,  it  is  not  the  variables  at  dummy  cells  that 
are  set  to  a  certain  value,  but  their  effect  on  the  implicit  operator  which  we  construct. 
Specifically,  the  boundary  conditions  are  folded  into  the  implicit  operator  by  making  a 
modification  to  the  Aitj  matrix  where  appropriate.  For  example,  the  inviscid  boundary 
conditions  at  the  wall  are  included  by  adding  the  matrix  Ci^Ew  to  Ait 2,  where  was 
defined  previously,  and  Ew  is  an  operator  that  reflects  the  momentum  normal  to  the  wall. 
Similarly,  the  approximate  implicit  viscous  Jacobian  Mn  is  modified  at  solid  boundaries. 
The  details  of  the  implicit  treatments  of  boundaries  are  described  in  Candler  (1988)  and 
Gokgen  (1989).  Note  that  the  viscous  wall  implicit  boundary  conditions  are  folded  into 
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the  block  matrix  Ai $  and  should  not  be  included  in  C 

6.  Other  Flux  Evaluation  Methods 

The  inviscid  flux  evaluation  method  discussed  earlier  is  one  of  many  methods  suitable 
for  gas  dynamic  simulations.  There  are  several  methods  used  widely  in  numerical  algo¬ 
rithms,  and  these  are  of  variable  degree  of  sophistication  and  complexity.  Different  flux 
evaluation  methods  exhibit  different  properties  in  terms  of  accuracy  and  computational 
cost.  Most  importantly,  the  amount  of  numerical  dissipation  that  is  associated  with  a  nu¬ 
merical  scheme  depends  largely  on  the  choice  of  flux  evaluation  and  solution  reconstruction 
methods.  In  general,  a  scheme  with  little  numerical  dissipation  produces  more  accurate 
results  on  a  given  computational  mesh. 

It  is  beyond  the  scope  of  this  work  to  make  comparisons  between  different  flux  eval¬ 
uations  methods.  However,  it  is  important  to  assess  how  well  the  numerical  method  of 
choice  performs  in  comparison  with  other  methods.  Therefore,  a  later  section  discusses 
in  more  detail  resuits  of  simulations  made  using  several  other  methods  of  evaluating  the 
inviscid  fluxes.  These  flux  evaluation  methods  are  not  covered  in  great  detail,  and  are  only 
presented  briefly. 

7.  Second  Order  Spatial  Accuracy 

Higher  order  accurate  inviscid  fluxes  can  be  obtained  through  appropriate  reconstruc¬ 
tion  of  the  solution  using  cell-centered  data.  Based  on  the  choice  of  the  reconstruction 
method,  the  accuracy  of  the  inviscid  fluxes  may  differ.  In  this  section,  we  briefly  present 
the  methods  used  in  this  work. 

The  flux  evaluation  method  discussed  previously  is  first  order  accurate  when  data 
from  adjacent  cells  are  used  to  calculate  the  flux.  The  method  can  become  second  order 
accurate  with  a  simple  modification.  Consider  the  Steger- Warming  flux  for  the  inviscid 
equations  in  one  spatial  dimension: 


F  =  F+(Ui)  +  F_(Ui+l)  . 

In  this  form,  the  split  fluxes  use  upwind  data  from  adjacent  cells.  A  higher  order  numerical 
flux  H  may  be  obtained  if  we  replace  the  arguments  of  F+  and  F_  with  a  left  state  U1' 
and  a  right  state  UR  as: 

H  =  F+{Ul)  +  F-{Ur) 

=  a+ul  +  a.ur  , 

where  UL  and  UR  are  extrapolated  conserved  quantities  based  on  pure  upwind  data.  It 
should  be  noted  that  this  simple  extrapolation  is  not  suitable  for  regions  of  the  flow  that 
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exhibits  large  gradients  in  U ,  such  as  shock  waves.  Therefore,  we  need  to  detect  the 
presence  of  a  shock,  and  use  no  extrapolation  in  that  case.  Our  implementation  uses 
pressure  jumps  across  cell  faces  in  order  to  detect  the  presence  of  shocks,  as  discussed 
earlier. 

The  MUSCL  variable  extrapolation  method  may  also  be  used  to  construct  higher 
order  Steger- Warming  and  modified  Steger- Warming  fluxes.  In  this  case,  the  left  and  right 
states  UL,  UR  are  constructed  based  on  non-linear  functions  of  the  ratio  of  adjacent  slopes 
of  U  computed  based  on  cell-centered  data.  Details  on  slope  limiters  and  their  use  in  the 
MUSCL  variable  extrapolation  can  be  found  in  Hirsch  (1991)  and  Laney  (1998). 
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IV.  Detailed  Vibrational  Relaxation  and  Dissociation  Modeling 

In  this  subsection  we  present  work  that  was  done  exclusively  on  nitrogen  flows.  This 
work  aims  to  address  modeling  issues  associated  with  vibrational  relaxation,  its  effect  on 
dissociation,  and  the  sensitivity  of  double-cone  flows  to  these  aspects  of  the  chemistry 
modeling.  We  do  not  concern  ourselves  with  the  numerical  treatment  of  the  equations 
that  are  solved,  and  the  focus  is  solely  on  the  physical  modeling. 

We  believe  that  from  a  modeling  perspective  it  is  easier  to  investigate  the  high  enthalpy 
nitrogen  cases.  Furthermore,  it  makes  sense  to  focus  our  efforts  on  the  nitrogen  flows,  for 
which  we  obtained  very  good  agreement  under  low  enthalpy  conditions.  Therefore,  in  the 
present  work  we  build  on  the  previous  modeling  efforts  to  study  the  CUBRC  LENS  high 
enthalpy  nitrogen  flows, 

We  show  results  of  vibrational  relaxation  in  nitrogen  using  the  latest  state-specific 
models  available  in  the  literature  for  vibrational  energy  relaxation  through  quantum  tran¬ 
sitions.  We  emphasize  that  although  the  equation  set  is  very  complex,  the  uncertainty  in 
the  rates  that  accompany  the  model  are  the  key  source  of  uncertainty  in  the  numerical 
predictions.  We  use  recently  developed  models  for  multi-quantum  state-specific  vibration- 
translation  (V-T)  and  vibration-vibration  (V-V)  relaxation  rates  in  simulations  of  the 
expanding  high  enthalpy  flow  in  the  nozzle.  Parametric  studies  of  the  sensitivity  of  the 
free-stream  conditions  to  modeling  of  the  V-T  and  V-V  relaxation  rates  are  presented.  The 
results  show  that  the  free-stream  conditions  are  nearly  insensitive  to  detailed  vibrational 
relaxation  modeling  under  these  conditions  and  that  with  possible  minor  modifications  the 
traditional  Landau- Teller  model  is  sufficient.  The  reason  for  the  discrepancies  observed 
between  the  present  high  enthalpy  experiments  and  simulations  is  not  explained  from  our 
findings  here. 

1.  Computational  Method 

The  standard  conservation  equations  have  been  modified  to  implement  the  state- 
specific  model  of  vibrational  relaxation.  In  the  case  of  state-specific  modeling  of  nitrogen 
molecules,  each  vibrational  state  is  treated  as  a  different  chemical  species,  and  thus  a 
separate  mass  conservation  equation  is  solved  for  each  vibrational  state.  The  diatomic 
species  are  ordered  in  terms  of  the  energy  that  is  associated  with  each  of  the  vibrational 
states.  For  nitrogen,  this  results  in  having  nv  +  1  mass  conservation  equations,  where  nv 
is  the  number  of  vibrational  quantum  states  considered;  the  additional  equation  is  due  to 
the  nitrogen  atoms  present  in  the  flow.  Note  that  we  do  not  include  all  4G  vibrational 
states  of  nitrogen,  because  the  equations  would  be  too  expensive  to  solve  on  an  adequate 
grid.  However,  using  fewer  levels  than  46  is  adequate  as  will  be  shown  in  subsequent 
subsections.  We  use  the  Park  rates  for  the  Arrhenius- type  dissociation  reactions,  and  the 
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reverse  reactions  are  calculated  using  Park’s  fits  for  the  equilibrium  constants.  The  details 
of  the  physical  modeling  are  discussed  in  detail  in  the  following  subsection. 

2.  Vibrational  State-Specific  Modeling 

One  of  the  major  uncertainties  in  the  prediction  of  high-enthalpy  nozzle  flows  is  how  to 
model  the  coupled  vibrational  relaxation  and  recombination.  As  the  gas  expands  through 
the  nozzle,  the  vibrational  modes  relax  by  vibration-translation  (V-T)  and  vibration- 
vibration  (V-V)  collisional  processes.  At  the  same  time,  the  atoms  recombine  to  form 
new  molecules.  These  molecules  may  be  formed  with  large  vibrational  energy;  this  en¬ 
ergy  is  then  redistributed  by  V-T  and  V-V  relaxation.  This  process  has  been  modeled  by 
many  different  authors,  including  Bray  (1968,  1970),  TVeanor  et  al.  (1968),  Center  and 
Caledonia  (1972),  Ruffin  (1995a,  1995b),  and  Josyula  and  Bailey  (2004),  This  work  shows 
that  anharmonic  effects  are  important,  as  well  as  the  competition  between  the  rates  of  the 
different  processes.  Ruffin’s  work  is  notable  in  that  full  master  equation  simulations  of 
hypersonic  nozzle  flows  show  excellent  agreement  with  experiments. 

In  this  work,  we  use  recently  developed  models  for  multi-quantum  state-specific  vibra¬ 
tion-translation  (V-T)  and  vibration-vibration  (V-V)  relaxation  rates  [Adamovich  (1997, 
1998)].  These  closed-form  rate  models  are  the  most  accurate  currently  available,  and  have 
been  compared  to  experimental  data  and  detailed  quantum  mechanical  calculations.  We 
model  the  expanding  gas  using  a  state-specific  approach,  representing  each  vibrational  state 
as  a  separate  chemical  species  at  each  point  in  the  flow.  Each  of  the  vibrational  species  can 
react  with  another  via  V-T  and  V-V  processes  as  the  gas  expands  in  the  nozzle.  We  also 
include  a  separate  species  for  nitrogen  atoms  to  allow  V-T  relaxation  of  the  N2  molecules 
by  atoms.  The  atoms  are  allowed  to  recombine  according  to  appropriate  rate  models. 

The  standard  nozzle  flow  analysis  code  uses  a  simple  Landau- Teller  relaxation  model 
for  the  vibrational  energy  modes.  Recombination  is  modeled  conventionally,  with  a  single 
rate  obtained  from  the  dissociation  rate  and  the  equilibrium  constant.  As  discussed  in 
previous  work,  we  use  the  Spalart-Allmaras  single-equation  turbulence  model  for  the  nozzle 
boundary  layer,  and  we  compute  from  the  reservoir  through  the  throat  to  the  test  section. 
For  the  LENS-Leg  I  nozzle  with  1.125"  inch  diameter  throat,  we  use  a  1007  x  128  grid 
which  has  been  shown  to  give  grid-converged  results. 

We  want  to  make  direct  comparisons  between  higher  fidelity  modeling  approaches 
and  the  standard  model.  Therefore,  we  use  the  same  grid  and  turbulence  model  for  all 
simulations.  The  forced  harmonic  oscillator  (FI IO)  model  of  Adamovich  et  al.  (1997)  is 
used  with  multi-quantum  V-T  and  V-V  transitions  allowed.  The  V-T  reaction  is  written 

N2{u  =  z)+M#N2(u  =  /)  +  M  (1) 
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where  M  is  any  collision  partner.  The  V-V  process  is 


n2(v  =  n)  +  n2(v  =  i2)  ^  n2(v  =  fx)  +  n2(v  =  f2)  (2) 


Thus  in  principle,  an  N2  molecule  can  transition  from  any  vibrational  state  to  any  other 
vibrational  state  by  either  of  these  mechanisms.  Therefore,  for  a  full  state-specific  model 
of  N2  (with  46  levels),  this  is  a  very  costly  simulation  on  an  adequate  grid.  We  can  reduce 
the  cost  of  the  simulation  by  only  considering  near-resonant  V-V  processes  and  limiting 
the  allowed  change  in  the  vibrational  quanta. 

The  modeling  of  the  recombination  process  is  unclear  because  we  must  specify  the 
nascent  vibrational  state  of  the  recombining  molecules.  Most  previous  work  assumes  that 
only  the  uppermost  vibrational  level  is  active  in  the  reactive  processes.  This  approach 
is  based  on  the  ladder-climbing  model  of  dissociation  which  is  not  relevant  in  hypersonic 
nozzle  flows. 


We  can  obtain  some  guidance  from  previous  quasi-classical  trajectory  analysis  [Bose, 
et  al.  (1997)]  of  the  exothermic  Zeldovich  reaction  02  +  N  — »  NO  -I-  O.  In  that  work, 
the  nascent  NO  vibrational  state  was  tracked  for  a  range  of  reactant  vibrational  and 
translational  states.  It  was  found  that  the  molecules  tend  to  be  formed  at  low  vibrational 
states.  To  a  first  estimate,  the  product  molecules  arc  equally  distributed  across  the  lower 
15  vibrational  quanta,  with  a  significant  reduction  in  the  population  beyond  that  level. 
Thus,  at  least  for  this  reaction,  the  molecules  are  more  likely  to  be  formed  at  the  low 
vibrational  levels.  With  this  observation  in  mind,  we  can  develop  two  different  modeling 
approaches  for  the  product  vibrational  state. 

First  consider  strict  detailed  balancing;  here  we  assume  that  the  equilibrium  constant 
for  each  state’s  population  must  hold.  In  this  case,  we  have  in  equilibrium  (denoted  by  *) 


K  t,A  =  Bffll  JML 

^  }  [N2(v)]-  [N2]*  [N2(t01* 


(3) 


Since  we  know  the  equilibrium  constant  for  the  N2  dissociation  and  the  equilibrium  Boltz¬ 
mann  distribution,  we  can  determine  /{^(v).  This  approach  effectively  distributes  the 
product  molecules  to  the  vibrational  energy  pool  based  on  the  local  Boltzmann  distribu¬ 
tion.  This  favors  population  of  the  low  levels,  and  it  clearly  satisfies  detailed  balance. 

A  second  approach  is  to  use  the  results  of  the  Zeldovich  reaction  analysis  discussed 
earlier,  and  assume  that  the  molecules  are  equally  disposed  to  the  vibrational  states.  Then 
the  expression  for  the  equilibrium  constant  becomes 


K^{v)  = 


[N]*[N]*  1 
[N2]'  nv 


(4) 
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where  nv  is  the  number  of  vibrational  states  being  modeled.  It  is  not  clear  that  such  a 
formulation  satisfies  detailed  balance;  however  if  the  V-T  and  V-V  rates  are  relatively  fast, 
the  gas  will  reach  an  equilibrium  governed  by  a  Boltzmann. 

2.  State-Specific  Modeling  Results 

Non-Reacting  Results 

First  let  us  consider  a  pure  N2  gas  expanding  from  the  reservoir  conditions  given 
by  T„  =  6720  K,  pa  =  8.680  kg/m3,  Cn  =  0.02647.  If  we  run  the  state-specific  code 
and  the  single  vibrational  temperature  Landau-Teller  model,  we  should  get  very  similar 
results.  As  shown  in  Table  1,  this  is  the  case  with  the  state-specific  code  giving  essentially 
the  same  results  as  the  baseline  code.  The  only  difference  should  be  due  to  using  the 
anharmonic  oscillator  model  in  the  state-specific  simulations.  In  addition,  the  vibrational 
energy  distribution  function  is  Boltzmann  at  the  nozzle  exit,  as  required  by  the  model. 
In  this  state-specific  simulation,  10  vibrational  levels  were  used;  we  will  show  that  this  is 
sufficient  to  represent  this  flow. 

Table  1  also  presents  the  results  for  state-specific  simulations  using  the  forced  harmonic 
oscillator  model  rate  expressions  for  single  quantum  V-T  transitions  and  for  a  full  multi¬ 
quantum  V-T  model  with  nearest  neighbor  V-V  transitions.  The  FHO  rates  are  faster  than 
the  Millikan  and  White  expression  for  the  Landau- Teller  relaxation  model,  so  the  state- 
specific  simulations  give  more  rapid  equilibration,  This  is  shown  in  Fig.  8,  which  plots  the 
FHO  and  Millikan  and  White  rates  for  the  u  =  ltou  =  0V-T  transition.  The  faster  FHO 
rates  result  in  a  lower  test-section  vibrational  temperature  and  a  correspondingly  higher 
velocity.  Somewhat  surprisingly,  the  multi-quantum  V-T  and  the  nearest- neighbor  V-V 
processes  do  not  have  an  appreciable  effect  on  the  vibrational  relaxation.  This  results  in 
essentially  no  difference  in  the  predicted  free-stream  conditions. 


Reacting  Results 

Initial  simulations  with  the  two  possible  recombination  models  showed  little  difference. 
Therefore,  the  detailed- balance  model  was  used  for  all  of  the  simulations  presented  here. 
Figure  8  shows  that  the  V-T  rate  for  N  atom  collisions  is  an  order  of  magnitude  faster 
than  the  molecule- molecule  rate.  Thus  we  expect  that  N  atoms  will  increase  the  vibrational 
relaxation  process.  Also,  as  the  atoms  recombine,  they  release  energy  which  increases  the 
flow  energy. 
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Table  1.  Simulation  results  for  non- reacting  nozzle  Hows  for  reservoir  conditions  corre¬ 
sponding  to  equilibrium  at  T0  =  6720  K,  pa  —  8.680  kg/m3  and  100%  N2.  Vibrational 
temperature  based  on  lowest  two  levels. 


Baseline 

State-Specific 

Landau-Teller 

St  ate -Specific 

FHO  Single-  Quantum  V-T 

State-Specific 

FHO  Multi- Quantum  +  V-V 

ti(m/s) 

3933 

3918 

3946 

3946 

p(e/™3) 

1.491 

1.503 

L481 

1.481 

T’(K) 

285.0 

296.7 

296.7 

296.7 

r„{K) 

3086 

3037 

2732 

2731 

PPitot 

23130 

23140 

23110 

23110 

1.0 

TO 

1.0 

1.0 

Table  2.  Simulation  results  for  reacting  nozzle  flows  for  reservoir  conditions  corresponding 
to  equilibrium  at  Ta  =  6720  K,  pa  —  8.680  kg/m3  and  97.35%  Na  by  mass. 


Baseline  State- Specific 

FHO  Multi-quant. 

State-Specific 

V-T  Multi-quant.,  16  levels 

State- Sped  fic 

V-T  +  V-V  [Aip=±n 

State-Specific 

V-T  +  V-V 

u(m/s) 

4108 

4137 

4149 

4127 

4126 

/j(g/m3)  1.407 

L388 

1.381 

1.393 

L393 

T{K) 

321.5 

344.6 

348.6 

342.6 

342.4 

r»(K) 

3168 

2535 

2544 

2536 

2536 

PPitot 

21900 

23840 

23850 

23800 

23800 

CN, 

0.9952 

0.9962 

0.9960 

0.9952 

0.9962 

Figure  8.  Plot  of  the  vibrational-translational  relaxation  rate  for  N2  (v  =  0)  — ♦  N2  ( v  — 
1)  computed  using  the  Landau- Teller  model  with  Millikan- White  fits  and  FHO  model. 
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Figure  9.  Distribution  of  the  vibrational  energy  among  states  for  a  10  level  calculation 
and  a  16  level  calculation  plotted  in  terms  of  the  mass  fraction  (left)  and  energy  fraction 
(right). 

As  shown  in  Table  2,  in  general,  the  reacting  results  are  similar  to  those  presented 
above.  However,  there  is  one  obvious  difference:  as  expected,  the  N  atoms  increase  the  re¬ 
laxation  rate  and  lower  the  predicted  test-section  vibrational  temperature.  The  additional 
energy  mostly  goes  into  increasing  the  free-stream  kinetic  energy. 

Interestingly,  the  V-V  rates  have  very  little  effect  on  the  flow  conditions.  Even  when 
V-V  transitions  that  have  Ai>  —  ±2,  there  is  little  effect  of  V-V  processes.  The  tabulated 
results  also  show  that  the  model  with  10  vibrational  levels  is  sufficient  to  capture  the  ener¬ 
getics  of  the  flow.  Increasing  the  number  of  levels  modeled  to  16  results  in  no  appreciable 
change.  This  result  is  further  quantified  in  Fig.  9,  which  plots  the  vibrational  population 
distribution  at  the  exit  of  the  nozzle  for  the  10  and  16  level  models.  Figure  9  also  plots 
the  energy- weighted  distribution  function  at  the  same  location;  note  that  almost  all  of  the 
energy  is  represented  by  the  first  few  levels  for  this  flow. 

These  results  are  somewhat  in  contrast  with  those  presented  in  previous  work.  How¬ 
ever,  they  are  the  first  to  use  the  recently-developed  forced  harmonic  oscillator  rates  which 
are  more  accurate  than  previous  rate  models.  Clearly  under  the  conditions  considered  here, 
the  nearest- neighbor  V-T  processes  dominate  the  recombination,  Also,  the  recombination 
of  atoms  without  favoring  the  upper  levels  has  a  significant  effect  on  the  vibrational  energy 
distribution  function  in  the  nozzle  flow. 

Interestingly,  Fig.  9  also  shows  that  the  vibrational  energy  modes  are  in  a  Boltzmann 
distribution  at  the  nozzle  exit.  This  is  contrary  to  much  of  the  previous  literature  on 
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hypersonic  nozzle  expansions  which  predict  a  complicated  vibrational  energy  distribution 
function.  (See  Ruffin  (1995a)  for  a  clear  discussion.)  Again,  this  is  a  result  of  the  dominant 
role  of  the  V-T  rates  predicted  by  the  FHO  rate  model  and  the  non- preferential  disposal 
of  the  recombination  energy. 

All  of  this  complicated  flow  modeling  results  in  a  relatively  simple  conclusion:  The 
simple  Landau-Teller  single  vibrational  temperature  model  will  work  well  for  high-enthalpy 
nozzle  flows  if  the  Millikan  and  White  relaxation  rates  are  adjusted  to  fit  the  modern  V- 
T  rate  data.  Of  course  this  conclusion  applies  to  pure  nitrogen  flows  only.  Ultimately, 
for  the  prediction  of  nozzle  conditions  in  a  hypersonic  wind  tunnel,  the  main  issue  is  the 
correct  prediction  of  the  rate  of  energy  release  from  the  vibrational  and  chemical  energy 
modes.  Subtle  differences  in  the  upper  vibrational  state  populations  have  little  effect  on 
the  macroscopic  flow  parameters  that  are  required  for  high-quality  aerodynamic  testing. 
Having  made  that  broad  generalization,  it  remains  to  be  seen  how  air  recombines  from  a 
high-enthalpy  reservoir  condition,  which  will  be  the  subject  of  future  studies. 

3.  Double-Cone  Flow  Results 

In  this  subsection  we  present  results  from  simulations  of  high  enthalpy  double-cone 
flows.  Following  the  discussion  of  the  detailed  nozzle  flow  simulations,  we  present  results 
obtained  under  several  of  the  computed  free-stream  conditions  presented  earlier,  as  well 
as  additional  cases  taken  from  Wadhams,  et  ai  (2003).  It  should  be  noted  that  all  of  the 
double-cone  predictions  presented  here  were  obtained  using  the  baseline  simple  harmonic 
oscillator  (SHO)  model  for  vibration  with  the  L-T  relaxation  model.  Table  3  summarizes 
the  free-stream  conditions  for  which  results  are  shown  in  this  subsection. 


Table  3.  Summary  of  free-stream  conditions  for  which  simulations  over  the  double-cone 
model  were  performed. 


Run  42 
(nominal) 

Run  42 
(computed) 
SHO  L-T 

Run  42 
(computed) 
FHO  v-t 

Run  42 
(computed) 
FHO  u-i  v-v 

Run  46 
{nominal) 

Run  46 
(computed) 
SHO  L-T 

Run  50 
(nominal) 

Run  40 
(nominal) 

3849 

4106 

4137 

4126 

3945 

4229 

3904 

3104 

P(g/m3) 

L468 

L407 

L388 

1.393 

1.958 

1,839 

1.509 

2.530 

T(K) 

268.7 

321.1 

344.6 

342.4 

281.7 

347.1 

272.5 

172.0 

T.(K) 

2947 

3170 

2535 

2536 

3072 

3083 

3143 

2617 

Cn, 

1. 0000 

0.9952 

0.9952 

0.9962 

0.9984 

0.9960 

0.9984 

1.0000 

The  initial  simulations  of  high  enthalpy  nitrogen  cases  (Runs  42  and  46)  presented 
in  Nompelis,  et  al.  (2003a)  were  performed  under  the  nominal  free-stream  conditions, 
but  also  under  free-stream  conditions  computed  with  CFD  by  expanding  the  flow  in  the 
LENS  nozzle  from  reservoir  conditions.  The  nominal  and  computed  conditions  were  very 
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similar,  and  simulations  under  both  sets  of  conditions  gave  very  similar  results  (less  than 
5%  difference  in  heat  transfer  rates  for  nitrogen).  For  the  purposes  of  this  discussion  we 
will  refer  to  these  results  as  “nominal.”  Recently,  the  reservoir  conditions  for  these  runs 
were  re-evaluated  as  presented  in  MacLean,  et  al.  {2005} ,  and  a  separate  set  of  reservoir 
conditions  was  proposed.  The  newly  proposed  free-strcam  conditions  suggest  an  increase 
of  nearly  10%  in  the  enthalpy  for  the  high  enthalpy  nitrogen  case  (Run  42).  All  results 
presented  here  were  obtained  from  the  proposed  reservoir  conditions,  where  the  free-strcam 
was  computed  via  nozzle  flow  analysis. 

Figure  10  plots  the  surface  pressure  and  heat  transfer  rate  to  the  model  for  Run  42. 
The  computed  conditions  from  Table  1  obtained  by  expanding  the  reservoir  conditions 
using  the  baseline  model  were  used  for  this  simulation.  Under  these  conditions,  the  size  of 
the  separation  zone  is  under-predicted  by  a  substantial  amount,  resulting  in  overall  poor 
agreement.  However,  agreement  is  better  for  the  pressure  in  the  attached  flow  region  along 
the  first  cone  and  inside  the  separation  zone. 

Figure  1 1  shows  heat  transfer  rate  results  computed  under  the  free-stream  conditions 
of  Table  2  obtained  using  the  more  physically  accurate  models:  the  state-specific  FHO 
model  with  multi-quantum  V-T  processes  but  neglecting  V-V  transitions,  and  the  state- 
specific  FHO  model  with  multi-quantum  V-T  relaxation  and  V-V  transitions  of  the  nearest 
±2  states.  The  baseline  result  is  also  plotted  for  reference.  The  results  worsen  when  the 
free-stream  conditions  obtained  using  the  more  accurate  model  are  used.  Also,  the  results 
are  mutually  consistent,  because  the  computed  free-stream  conditions  of  Table  2  are  very 
similar. 


Figure  10.  Numerical  predictions  and  experimental  measurements  for  Run  42.  The 
computed  conditions  using  the  baseline  model  from  Table  2  were  used  for  this  simulation; 
the  results  under  nominal  conditions  are  plotted  for  reference. 
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Figure  11.  Comparison  of  predicted  heat  transfer  rate  for  Run  42  using  the  state-specific 
FHO  multi-quantum  V-T  model,  and  the  state-specific  FIIO  multi-quantum  V-T  model 
with  V-V  transitions  for  the  nearest  ±2  states.  Results  under  nominal  conditions  are 
plotted  for  reference. 

Figure  12  shows  results  for  Run  46  conditions  computed  from  the  new  reservoir  con¬ 
ditions.  The  size  of  the  separation  zone  is  substantially  under- predicted  under  the  new 
conditions  obtained  with  the  baseline  model.  However,  the  level  of  pressure  along  the  first 
cone  and  in  the  separated  region  matches  the  measured  values.  There  is  an  overall  increase 
in  the  heat  transfer  rate  predictions,  which  is  closer  to  the  experimental  values  along  the 
first  cone,  but  the  agreement  is  still  poor  at  reattachment.  It  is  interesting  to  note  that 
although  the  predicted  heat  transfer  rate  increased  consistently  along  the  entire  surface  of 
the  model,  there  is  still  a  large  difference  at  the  peak  and  after  reattachment.  Also,  if  we 
examine  the  shear  layer  that  forms  at  the  edge  of  separation  and  characterize  it  in  terms 
of  a  Reynolds  number  as  defined  in  Birch  and  Keyes,  we  see  that  for  Run  46  this  quantity 
runs  dangerously  close  to  the  transition  criterion  specified  by  the  authors.  They  compute 
the  free  shear-layer  Reynolds  number  based  on  the  high-speed  side  as:  Rca  —  pvs/n,  where 
s  is  the  length  along  the  shear  layer.  This  is  shown  in  Fig.  13,  which  plots  the  shear  layer 
Reynolds  number  for  this  case  along  with  Runs  42,  43  and  the  low-enthalpy  case  Run  35. 
Although  this  dimensionless  quantity  is  below  the  transition  criterion,  it  has  a  trend  the 
indicates  the  flow  might  be  unsteady.  Based  on  these  results,  we  believe  that  it  is  likely 
that  this  case  is  transitional,  and  thus  it  is  not  useful  for  the  purposes  of  this  study 
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Figure  12.  Numerical  predictions  and  experimental  measurements  for  Run  46.  The 
computed  conditions  using  the  baseline  model  from  Table  2  were  used  for  this  simulation; 
the  results  under  nominal  conditions  are  plotted  for  reference. 


Figure  13.  Free  shear-layer  Reynolds  number  for  Run  46,  Run  42,  Run  43  (air  case)  and 
Run  35  (low  enthalpy  case). 
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Figure  14.  Computed  surface  pressure  and  heat  transfer  rate  for  Run  50  (blunt  double¬ 
cone  in  10  MJ/kg  nitrogen). 

We  also  simulated  the  flow  over  a  blunted  double-cone  model  in  a  high  enthalpy 
nitrogen  flow  (Run  50).  Run  50  is  nominally  the  same  as  Run  42,  and  was  included  in  the 
LENS  experiments  to  assess  the  effect  that  a  blunted  nose  has  on  the  flow.  The  blunted 
nose  has  a  strong  stagnation  at  the  model  tip,  which  forms  a  large  entropy  layer  that 
extends  into  the  separation  zone.  Chemical  effects  in  the  entropy  layer  may  possibly  alter 
the  flowfield,  therefore  it  is  important  to  know  how  this  effect  changes  the  basic  flowfield. 
The  simulation  results  along  with  the  experimental  data  for  this  case  are  plotted  in  Fig.  14. 
We  see  good  agreement  for  the  surface  pressure  and  heat  transfer  rate  along  the  first  cone, 
however  the  size  of  the  separation  is  under-predicted  substantially.  This  is  a  surprising 
result  because  the  free-stream  and  reservoir  conditions  for  Run  50  are  nearly  identical  to 
Run  42  and  the  CFD  results  for  both  cases  are  very  similar.  However,  this  is  not  reflected 
in  the  experimental  measurements. 

Finally,  we  simulated  the  flow  over  the  double-cone  model  under  the  nominal  free- 
stream  conditions  of  Run  40,  which  is  a  moderate  enthalpy  case  nominally  at  5  MJ/kg. 
This  case  exhibits  large  scale  unsteadiness  as  predicted  by  the  simulation,  in  contrast  to 
what  is  seen  in  the  experiments.  This  is  shown  in  Fig.  15,  which  plots  the  heat  transfer 
rate  to  the  model  for  Run  40.  The  separation  zone  extends  almost  over  the  entire  length 
of  the  geometry  and  is  unsteady. 
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Figure  15.  Computed  surface  heat  transfer  rate  for  Run  40. 

We  have  simulated  high  enthalpy  experiments  in  nitrogen  performed  in  the  LENS 
facility  using  the  LENS  Leg  1  tunnel.  We  have  studied  the  experiments  by  performing 
simulations  of  the  nozzle  flow,  expanding  the  reservoir  conditions  from  the  throat  into  the 
test-section.  We  studied  Run  42  extensively  using  the  most  recently  developed  models  for 
vibrational  relaxation,  including  V-T  and  V-V  relaxations  processes.  We  demonstrated 
that  for  the  Hows  of  interest  the  free-stream  conditions  are  not  very  sensitive  to  the  mod¬ 
els  of  the  relaxation  process  and  that  the  standard  SIIO  model  with  the  L-T  model  for 
relaxation  gives  accurate  results. 

We  also  investigated  the  effect  of  bluntness  on  the  high  enthalpy  double-cone  flow 
by  simulating  two  different  experiments  at  the  same  conditions  for  a  sharp  and  a  blunted 
double-cone  model.  The  CFD  results  were  very  similar  for  the  two  experiments,  but  this 
was  not  reflected  in  the  experimental  measurements.  The  size  of  the  separation  zone  was 
larger  for  the  blunted  case.  The  reason  for  this  difference  is  unknown,  and  a  more  detailed 
analysis  is  needed  to  answer  this  question,  as  well  as  additional  experiments  to  verify  this 
trend  experimentally. 

Finally,  from  this  study  we  were  unable  to  identify  the  main  reason  for  the  discrep¬ 
ancies  between  experiments  and  simulations  that  was  observed  initially,  and  in  fact  with 
the  proposed  reservoir  conditions  the  agreement  was  worse.  We  believe  that  the  proposed 
reservoir  conditions  are  not  accurate. 


42 


V,  Nozzle  Flow  Analysis  with  Non-intrusive  Diagnostics 

In  view  of  the  poor  agreement  that  we  have  observed  for  high  enthalpy  air  experi¬ 
ments  with  double-cone  geometries,  we  shifted  focus  to  analyzing  the  nozzle  flow  for  strong 
nonequilibrium  effects  that  may  result  in  incorrectly  inferring  the  free-stream  conditions. 
The  experimental  group  has  invested  a  lot  of  resources  in  establishing  a  non-intrusive 
methodology  for  making  measurements  inside  the  nozzle  for  certain  flow  quantities,  such 
as  nitric  oxide  (NO)  content.  In  this  subsection  we  present  work  that  was  done  in  an 
effort  to  interpret  the  recent  NO  content  measurements  made  at  LENS.  It  should  be  noted 
that  the  levels  of  NO  that  are  reported  by  the  experimental  group  in  Parker  et  al.  (2006) 
are  lower  than  we  predict  with  CFD  calculations  under  the  same  conditions.  This  work 
summarizes  the  methodology  and  key  findings. 

In  this  work  we  re-examine  the  nozzle  flow,  including  the  newest  runs,  to  understand 
these  measurements.  To  do  so,  we  study  the  chemical  reaction  mechanisms  that  take 
place  in  the  nozzle  in  order  to  explain  the  reason  behind  the  elevated  level  of  NO  that 
is  predicted  with  CFD  compared  to  what  is  measured  experimentally.  We  compare  the 
physical  models  that  we  employ  to  the  more  detailed  models  employed  by  Park  for  conical 
nozzle  flow  simulations  at  similar  conditions  [Park  (2006)]. 


1.  Modeling  Assumptions  and  Overview 

The  high  enthalpy  LENS  experiments  performed  in  support  of  the  code  validation 
study  were  done  at  nominally  5,  10  and  15MJ/kg  free-stream  specific  enthalpies.  These 
conditions  correspond  to  flight  Mach  numbers  of  about  10,  15,  and  17  respectively.  At 
these  conditions  (above  5MJ/kg),  when  the  air  gas  stagnates,  the  temperature  rises  to 
high  levels,  where  dissociation  and  internal  energy  excitation  take  place.  These  processes 
are  modeled  in  CFD.  The  models  employed  in  CFD  codes  and  for  inferring  the  free-stream 
conditions  make  use  of  published  data.  Models  are  employed  for  calculating  the  rates  at 
which  energy  relaxation  among  different  modes  takes  place.  Models  are  also  employed  to 
account  for  the  distribution  of  energy  in  the  different  modes  for  each  species,  and  how 
the  different  modes  interact  with  each  other.  Although  the  physical  models  may  allow 
for  different  energy  modes  to  be  out  of  equilibrium,  the  basic  assumption  is  that  all  the 
states  for  a  given  energy  mode  are  in  equilibrium  and  that  a  Boltzmann  distribution 
exists.  Vibrational,  electronic  and  rotational  are  the  energy  modes  that  are  relevant  for 
the  conditions  of  interest  in  addition  to  the  translational  mode. 

In  this  section  we  discuss  the  modeling  assumptions  that  support  the  methodology 
employed  in  CFD,  and  the  calculations  of  the  nozzle  flow.  Vibrational  relaxation,  chem¬ 
ical  reaction  rates,  dissociation,  recombination  and  the  distribution  of  nascent  states  for 
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recombining  molecules  are  discussed. 


Tkermochemical  Equilibrium  in  the  Reservoir 

In  the  LENS  experiments,  shock  waves  heat  the  test  gas  to  high  enthalpy,  and  then 
the  gas  expands  through  the  throat  into  the  nozzle  to  high  Mach  number.  During  this 
process,  the  gas  spends  sufficient  time  in  the  reservoir  to  reach  thermochemical  equilibrium. 
When  the  stagnant  gas  reaches  equilibrium  at  high  enthalpy  conditions,  part  of  the  total 
enthalpy  is  stored  in  the  form  of  chemical  energy  via  the  existence  of  radicals  (atomic 
species)  and  nitric  oxide  (NO).  Other  species  -  such  as  N2O  -  are  also  created  but  they 
exist  in  negligible  amounts  under  these  conditions,  Internal  energy  modes  of  the  gas  are 
excited,  and  a  portion  of  the  energy  is  stored  in  vibrational  and  electronic  energy  modes. 
With  increasing  total  enthalpy,  a  larger  portion  of  the  energy  is  stored  in  chemical  form 
when  the  gas  is  stagnant  and  in  thermochemical  equilibrium.  Table  4  shows  how  the 
total  energy  is  distributed  among  the  different  energy  modes  when  the  gas  is  stagnant  at 
different  conditions. 


Table  4.  Distribution  of  the  total  enthalpy  in  the  reservoir  for  different  conditions  at 
p  =  500  atm. 


Hq 

5  MJ/kg 

10  MJ/kg 

15  MJ/kg 

T(  K) 

3780 

6085 

7897 

Thermal 

77.6% 

64.7% 

58.3% 

Vibrational 

13.8% 

11.8% 

9.6% 

Chemical 

8.6% 

23.4% 

32.1% 

The  test-gas  in  the  reservoir  remains  in  equilibrium  at  a  state  that  is  assumed  not 
to  change  with  time  throughout  the  useful  test-time  of  the  experiment.  Measurements 
are  made  in  the  reservoir  to  help  measure  the  useful  test-time  window.  The  shock  heated 
gas  expands  to  high  Mach  number  in  the  nozzle.  During  the  sudden  expansion,  the  gas 
is  in  a  nonequilibrium  state  and  energy  relaxation  and  recombination  of  the  dissociated 
species  occur  in  the  nozzle.  This  process  is  not  well  understood  and  we  believe  that  it  is 
not  adequately  described  by  the  air-chemistry  models  that  are  presently  employed.  This 
is  the  subject  of  investigation  in  the  present  work.  The  chemical  reaction  rates  and  their 
underlying  assumptions  are  discussed  below. 


Vibrational  Energy  Modes  and  Relaxation 
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It  has  been  shown  that  vibrational  energy  relaxation  plays  an  important  role  for  the 
low  enthalpy  LENS  experiments  [Nompelis  (2003)].  Vibrational  energy  is  accounted  in 
the  CFD  methodology  by  including  a  vibrational  energy  equation  that  is  solved  at  every 
point  in  the  flow.  The  distribution  of  vibrational  energy  is  described  by  a  simple  harmonic 
oscillator  (SHO)  model  for  each  of  the  diatomic  species.  In  reality,  the  vibrational  mode 
has  a  number  of  quantum  states  that  are  characteristic  for  each  species.  A  conservation 
equation  may  be  considered  for  the  population  of  each  state,  but  such  a  calculation  would 
required  accurate  knowledge  of  the  state-to-state  and  state-specific  V-T  transitions  and 
would  result  in  a  large  set  of  conservation  equations.  Because  most  of  the  energy  is 
concentrated  on  the  first  excited  state,  the  SHO  model  appears  to  be  adequate  for  these 
flows.  Detailed  studies  of  vibrational  relaxation  for  high  enthalpy  nitrogen  flows  in  the 
LENS  facility  accounting  for  neighboring  and  near- neighboring  state  transitions  have  been 
performed.  These  simulations  showed  that  the  detailed  models  had  little  effect  on  the 
computed  free-stream  conditions  [Nompelis  (2005)].  For  vibrational  relaxation  with  the 
translational-rotational  mode  (V-T)  we  assume  the  Landau- Teller  model  is  valid. 

Chemical  Reactions 

Chemical  reactions  take  place  when  the  flow  is  in  chemical  nonequilibrium.  The 
chemical  reactions  that  are  important  in  these  flows  are  dissociation  and  recombination  of 
diatomic  species  (N2,  O2,  NO)  which  take  the  form: 

Nj  +  MfiN  +  N  +  M 

O2  4*  M  -  O  +  O  4-  M 

NO  4  M  ^  N  4  O  4  M 

where  M  is  a  generic  collision  partner,  and  the  Zeldovich  reactions: 

N2  4-  O  NO  4-  N 
NO  4-  0  02  4-  N 

The  total  rate  of  reaction  for  each  of  these  processes  is  the  sum  of  the  forward  and  the 
backward  reaction  rates.  The  reaction  rate  coefficients  k j  and  kb  that  appear  in  the  rate 
expressions  are  functions  of  temperature,  they  are  different  for  each  reaction,  and  in  the 
case  of  dissociation  reactions  they  are  different  for  each  collision  partner  M.  Expressions 
for  the  rate  coefficients  are  taken  from  the  literature.  The  parameters  that  appear  in 
the  functional  form  of  the  rate  constants  are  typically  obtained  experimentally  and  are 
valid  for  a  range  of  temperatures.  Typically,  only  the  form  of  the  forward  reaction  rate 
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coefficient  kj  is  given  in  the  literature.  The  backward  rate  coefficient  is  calculated  by  using 
the  principle  of  detailed  balance. 


It  is  assumed  that  the  reaction  moves  in  one  direction  or  another  based  on  what 
the  composition  and  the  state  of  the  gas  are  relative  to  equilibrium.  The  composition 
at  equilibrium  is  a  thermodynamic  property  of  a  mixture  of  pure  products  and  reactants 
involved  in  a  reaction  and  is  a  function  of  the  thermodynamic  state  of  the  mixture.  Using 
the  equilibrium  constant  Keq  the  following  relationship  is  used  as  a  result  of  detailed 
balance: 

K  —  — 

Keq  ~  kb 


Published  data  exist  for  determining  the  equilibrium  constant  It  should  be  noted 

that  the  principle  of  detailed  balance  holds  under  the  assumption  that  the  reaction  is  a 
one-step  process  from  products  to  reactants.  Therefore,  the  calculation  of  the  backward 
rate  coefficient  is  only  phenomenological. 


Coupling  of  Reactions  and  Energy  Modes 

Because  the  chemical  reactions  are  a  result  of  collisions  between  particles,  the  reaction 
rate  coefficients  are  functions  of  temperature.  However,  in  the  case  of  dissociation,  a  given 
dissociating  molecule  may  be  in  any  vibrational  state.  Thus,  if  diatomic  molecules  are  on 
average  at  a  high  vibrational  energy  (elevated  Tv),  on  average,  more  collisions  will  result 
in  dissociation,  and  this  is  because  more  molecules  exist  close  to  their  dissociation  limit. 
Several  models  have  been  formulated  to  account  for  this  effect.  In  the  present  work  the 
model  of  Park  is  used  which  modifies  the  temperature  at  which  the  forward  reaction  rate 
coefficients  are  evaluated  to  use  a  geometric  average  \/TTv.  This  is  a  simple  means  of 
accounting  for  the  effect  of  elevated  vibrational  energy  in  dissociation.  The  backward  rate 
coefficient  for  dissociation  reactions  is  not  calculated  based  on  the  geometric  temperature 
and  thus  it  does  not  depend  on  Tv. 

When  atoms  recombine,  a  molecule  is  created  with  a  finite  amount  of  energy.  The 
energy  of  the  particle  accounts  for  the  sum  of  the  translational-rotational,  electronic  and 
vibrational  energies.  Because  we  use  a  macroscopic  description  of  the  recombination  pro¬ 
cess,  we  assume  that  the  nascent  molecules  are  created  at  all  possible  vibrational  states 
following  an  equilibrium  distribution  at  the  local  Tv.  Since  electronic  energy  is  neglected, 
the  remaining  energy  is  assumed  to  be  distributed  in  translation-rotation  in  an  equilibrium 
distribution. 

It  has  been  speculated  that  molecules  may  recombine  at  excited  electronic  states.  In 
particular,  long-lived  diatomic  and  atomic  oxygen  at  excited  electronic  states  may  exist 
in  the  reservoir.  These  energetic  particles  may  interact  with  other  species  such  that  they 
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participate  in  recombination  and  the  Zeldovich-type  reactions  in  a  preferential  direction. 
These  effects  have  not  been  considered  in  previous  studies  of  the  LENS  experiments. 


2.  CUBRC  LENS-I  Nozzle  Flow  Diagnostics 

Parker  et  al.  used  infrared  laser  energy  absorption  to  measure  the  velocity,  temperature 
and  nitric  oxide  mass  fraction  at  the  exit  plane  of  the  CUBRC  LENS-I  nozzle  [Parker  et  al, 
(2006),  (2007) j  The  NO  rotational  transition  at  5447nm  was  used  for  the  measurements. 
The  Doppler  shift  of  the  absorption  feature  gives  an  accurate  means  of  measuring  the  gas 
velocity,  while  the  width  of  the  feature  gives  an  estimate  of  the  gas  temperature  using  the 
HITRAN  2K  database  of  line  shapes.  Then  with  the  temperature  known,  the  absolute 
concentration  of  nitric  oxide  can  be  computed  with  the  assumption  that  the  rotational 
and  vibrational  modes  are  in  a  Boltzmann  distribution  at  the  measured  temperature.  We 
will  focus  on  Parker’s  h0  —  10.5MJ/kg,  p0  =  54.56  MPa  test  case,  for  which  he  measured 
an  exit  plane  velocity  of  4517  m/s,  temperature  in  the  range  of  250  —  300  K,  and  NO  mass 
fraction  of  1.5%.  A  quick  analysis  shows  that  these  measurements  imply  a  very  small  level 
of  internal  energy  storage  in  the  test-section  gas.  The  total  enthalpy  of  the  gas  mixture  is: 


ho  —  CpT  +  ev  “i-  echeTn  ^  u 
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Here  l/2u2  =  10.20  MJ/kg,  and  with  cp  ~  lOQGJ/kgK  and  T  =  275  K,  we  have  just 
25kJ/kg  remaining  in  the  internal  energy  modes.  With  1.5%  NO  mass  fraction  (h^0  — 
3.0 MJ/kg),  the  remaining  vibrational  and  chemical  energy  is  essentially  zero.  Therefore, 
according  to  Parker’s  measurements,  the  test-section  gas  has  a  negligible  level  of  vibrational 
and  chemical  energy. 

Presumably  the  differences  in  the  double-cone  measurements  and  simulations  are  due 
to  problems  with  the  predicted  conditions  at  the  exit  of  the  CUBRC  LENS-I  nozzle. 
However,  MacLean  et  al.  (2007)  show  some  improvement  in  the  double-cone  results  with 
the  use  of  the  CVDV  model  for  vibration-dissociation  coupling.  In  any  case,  it  is  important 
to  try  to  rationalize  Parker’s  results  with  nozzle  flow  field  models.  Therefore,  the  goal  of 
this  work  was  to  reassess  the  nozzle  flow  and  thermo-chemical  modeling  approaches  to  try 
to  improve  the  prediction  of  the  free-stream  conditions. 

Recently,  Park  did  a  complex,  vibrational  state-to-state  analysis  of  high-enthalpy  noz¬ 
zle  flows  [Park  (2006)].  His  approach  solves  the  vibrational  master  equation,  including 
state-specific  dissociation  and  recombination  processes.  The  rates  are  computed  with  the 
most  accurate  approaches  available.  Park’s  work  is  currently  the  most  advanced  high- 
enthalpy  nozzle  analysis.  It  showed  a  number  of  important  things  about  these  flows.  First, 
vibrational  non-equilibrium  is  not  an  important  factor  in  air.  The  presence  of  a  large 
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number  of  oxygen  atoms  in  the  gas  causes  very  rapid  relaxation  of  the  vibrational  modes. 
(It  is  critically  important  to  use  the  measured  data  for  the  N2-O  and  O2-O  relaxation 
times,  rather  than  the  standard  Millikan  and  White  relaxation  times.)  Also,  according  to 
Park’s  analysis  there  is  not  a  significant  overpopulation  of  the  upper  vibrational  states. 
Similarly,  the  vibrational  energy  modes  of  the  diatomic  species  are  quite  close  to  one  an¬ 
other,  with  O2  relaxing  slightly  more  slowly  than  N2.  However,  from  an  energy  budget 
perspective,  this  difference  is  inconsequential.  (Remember  that  our  goal  is  to  track  the 
overall  energy  budget  of  these  flows  -  thus  full  details  are  not  required  unless  they  affect  the 
energy  relaxation  and  redistribution  process.)  Finally,  Park’s  results  for  ha  —  lOMJ/kg, 
p0  —  1000  atm  (101.3  MPa)  case  with  similar  A/ A*  as  Parker’s  LENS-I  case  gives  a  nitric 
oxide  mass  fraction  of  3.92%  -  significantly  higher  than  measured  by  Parker,  Park’s  higher 
reservoir  pressure  would  presumably  increase  the  recombination  rates  relative  to  that  of 
Parker,  and  should  therefore  be  lower  than  Parker’s  1 .5%. 


3.  Nozzle  Flow  Model 

We  use  a  conventional  single  vibrational  energy,  thermo-chemical  nonequilibrium  ki¬ 
netics  model  for  the  nozzle  flow.  We  use  5  chemical  species  (N2,  O2,  NO,  N  and  O) 
and  neglect  the  effect  of  ions.  The  kinetics  model  is  from  Park  (1988)  and  we  use  the 
same  set  of  vibrational  relaxation  rate  correction  factors  given  by  Park  to  match  the  avail¬ 
able  experimental  data.  We  solve  the  axisymmetric  Navier-Stokes  equations  along  with 
the  Spalart-Allmaras  turbulence  model  with  a  compressibility  correction.  More  details 
about  the  nozzle  simulation  approach  are  available  in  Ref.  14,  The  nozzle  code  has  been 
used  to  predict  a  wide  variety  of  CUBRC  reflected  shock  tunnel  nozzle  flows  [MacLean  et 
al.  (2005)]. 

This  thermo-chemical  model  cannot  capture  vibrational  energy  inversions,  disparate 
vibrational  energy  distributions  between  the  diatoms,  and  the  details  of  vibration-dissoc¬ 
iation  /  recombination  interactions.  Therefore,  we  must  first  validate  this  approach  to 
show  that  the  model  captures  the  key  energy  redistribution  processes  during  the  expansion 
of  high-enthalpy  air.  We  compare  to  the  lOMJ/kg,  lOOOatm  case  of  Park;  we  use  a 
conical  nozzle  of  5  meter  length  and  an  exit-to- throat  area  ratio,  A/A* ,  of  1000.  For  this 
comparison,  we  treat  the  flow  as  inviscid  so  that  the  area  variation  of  the  nozzle  is  not 
affected  by  the  boundary  layer  displacement.  In  addition  to  the  full  CFD  simulation,  wc 
compare  to  a  quasi-one- dimensional  nozzle  calculation.  This  method  includes  the  complete 
thermo-chemical  model  used  in  the  CFD,  and  simply  integrates  the  governing  equations 
from  the  throat  to  the  exit  plane  of  the  nozzle.  Figure  16  compares  the  CFD  and  quasi-one- 
dimensional  simulation  results  for  the  Park  case.  Note  the  excellent  agreement  between  the 
two  methods.  Interestingly,  the  non- uniformities  in  full  CFD  simulation  cause  variations 
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in  the  temperature  and  chemical  state,  however  the  final  exit  conditions  are  very  similar. 

The  nozzle  exit  conditions  are  compared  with  those  of  Park  in  Table  5;  note  the  ex¬ 
cellent  agreement  with  the  main  flow  parameters.  The  present  simulations  predict  slightly 
lower  vibrational  temperatures  than  Park’s  detailed  simulation.  Also,  we  predict  a  larger 
NO  mass  fraction  than  Park.  This  is  not  altogether  surprising  since  there  are  substantial 
differences  in  the  model  and  rates  used  by  Park.  Thus,  given  the  large  differences  in  the 
modeling  approach,  this  agreement  is  very  reasonable.  They  are  also  consistent  with  the 
comparisons  with  one-  and  two-temperature  nozzle  simulations  shown  by  Park. 


Table  5.  Computed  conditions  for  Park’s  lOMJ/kg  nozzle  flow  case. 


Park 

CFD 

Quasi-l-D 

T(K) 

674 

684 

744 

TA  K) 

990 

920 

936 

u  (m/s) 

4243 

4260 

4246 

P  (g/™3) 

11.17 

11.11 

11.02 

CNO 

3.93 

4.86 

4.86 

Co2 

20.26 

20.34 

20.65 

Table  6.  Computed  conditions  for  Parker’s  10.5MJ/kg  case  (Run  271). 


Parker 

CFD 

r(K) 

250-300 

620 

TA  K) 

930 

u  (m/s) 

4517 

4354 

CNO 

1.5 

5.42 

COj 

20.26 

19.10 

Co 

1.30 
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Figure  16.  Comparison  of  full  CFD  and  quasi-one-dimensional  simulations  of  the  center¬ 
line  conditions  for  Park’s  lOMJ/kg  test  case.  Right  figure  shows  an  enlargement  near  the 
throat. 


4.  Baseline  CUBRC  LENS-I  Nozzle  Simulation 

We  used  the  standard  thermo-chemical  model  to  simulate  the  LENS-I  nozzle  flow  at 
conditions  corresponding  to  the  10.5MJ/kg  and  56.54  MPa  run  of  Parker  (Run  271).  The 
resulting  nozzle  centerline  flow  field  is  shown  in  Figure  17,  and  the  exit  plane  conditions 
are  compared  with  the  measurements  of  Parker  in  Table  6.  Note  that  there  is  little  vibra¬ 
tional  nonequilibrium  in  this  case,  and  that  the  gas  composition  freezes  near  the  nozzle 
throat.  We  predict  an  NO  mass  fraction  of  5.42%,  significantly  above  the  1.5%  of  Parker. 
Correspondingly,  the  CFD  gives  a  lower  velocity  and  a  higher  temperature  than  measured. 
Thus,  the  CFD  predicts  that  a  significant  amount  of  energy  is  stored  in  the  vibrational  and 
chemical  energy  modes  (33MJ/kg  in  vibration  and  363kJ/kg  chemical  energy,  compared 
to  approximately  45MJ/kg  in  chemical  energy  for  Parker’s  data). 

This  is  a  very  significant  discrepancy,  and  is  certainly  a  potential  reason  for  the  poor 
comparison  with  the  double-cone  data  for  high  enthalpy  air.  In  the  remainder  of  the  section 
we  investigate  potential  reasons  for  this  difference.  It  should  be  noted  that  our  model  does 
not  include  the  full  state- to-state  kinetics  as  in  Park.  However,  based  on  the  results  of 
Park  and  the  comparison  presented  above,  it  is  unlikely  that  such  a  model  would  shed  light 
on  this  discrepancy.  Park  did  not  show  any  appreciable  vibrational  population  inversion 
or  other  effect  that  could  invalidate  the  present  modeling  approach.  Also,  because  of  the 
rapid  vibrational  relaxation  in  the  nozzle,  it  is  unlikely  that  vibration  plays  an  important 
role  in  the  chemical  kinetics.  Rather,  we  focus  on  additional  issues  that  are  not  considered 
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in  conventional  nozzle  flow  models. 


Figure  17.  Computed  centerline  flow  conditions  for  Parker’s  10,5MJ/kg  case  using  the 
baseline  thermo-chemical  model. 

5.  Nitric  Oxide  Destruction  Process 

The  main  issue  raised  by  the  comparison  with  Parker’s  data  is  the  excessive  level  of 
nitric  oxide  at  the  nozzle  exit  plane.  Let  us  consider  how  NO  is  consumed  as  the  gas 
expands  from  the  reservoir  (where  the  NO  mass  fraction  is  13.05%).  There  are  three 
reactions  involving  NO.  The  dissociat ion/ recombination  of  NO: 

N0  +  M^N  +  0  +  M 

and  the  Zeldovich  reactions: 

N2  +  O  ^  NO  +  N 

NO  +  Ov^Oj+N 

These  reactions  are  written  so  that  they  are  endothermic  as  they  proceed  to  the  right. 
During  the  expansion  of  the  flow  in  the  nozzle,  the  first  reaction  proceeds  to  the  left,  pro¬ 
ducing  additional  NO.  The  first  Zeldovich  reaction  consumes  NO  and  N  atoms  (proceeding 
to  the  left)  in  a  rapid  exothermic  reaction.  The  second  Zeldovich  reaction  also  destroys 
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NO  in  an  endothermic  reaction  that  produces  N  atoms,  which  are  then  used  in  the  first 
Zeidovich  reaction.  Therefore,  the  rate-limiting  step  is  the  availability  of  N  atoms.  Initially 
there  are  N  atoms  in  the  reservoir,  but  these  are  rapidly  consumed  by  the  recombination 
reaction  and  the  first  Zeidovich  reaction.  Once  they  are  used  up,  the  first  reaction  cannot 
take  place  without  N  atoms  being  formed  by  the  second  Zeidovich  reaction. 

Figure  18  illustrates  this  process  more  quantitatively  by  plotting  the  rate  of  NO  for¬ 
mation  as  a  function  of  distance  for  each  of  the  above  reactions.  We  plot  the  quantity 
ujno fpu,  where  iuno  is  the  chemical  source  term  for  nitric  oxide  due  to  each  of  the  three 
reactions.  This  gives  the  non-dimensional  rate  of  NO  formation  per  unit  length  of  the 
nozzle.  In  this  LENS-I  nozzle,  the  sonic  throat  is  located  at  approximately  0.154  m  from 
the  x  =  0  location;  at  this  point,  the  temperature  begins  to  change  rapidly,  and  the  gas  be¬ 
gins  its  recombination  process.  Note  that  initially  the  dissociation/ recombination  reaction 
produces  NO,  while  the  two  Zeidovich  reactions  consume  the  newly  formed  NO.  Then,  as 
the  gas  continues  to  cool  the  Zeidovich  reactions  quickly  consume  the  NO,  while  the  re¬ 
combination  reaction  produces  only  a  small  amount.  Note  that  the  two  Zeidovich  reactions 
produce  NO  at  essentially  the  same  rate;  this  confirms  that  the  second  reaction  is  the  rate 
limiting  step  (especially  late  in  the  expansion  process).  Interestingly,  at  x  =  0.18  m  there 
is  a  sudden  reversal  of  the  NO  destruction  process  due  to  a  recompression  that  increases 
the  centerline  temperature.  The  Zeidovich  reactions  then  continue  to  consume  NO,  and 
by  about  x  =  0.4  m  the  gas  has  cooled  and  accelerated  so  much  that  the  reactions  arc 
essentially  frozen.  Very  little  reaction  occurs  in  the  remainder  of  the  nozzle. 

Figure  18  shows  that  there  is  a  complicated  balance  between  the  NO  reactions,  the 
rates  of  formation  of  the  reactants  in  these  reactions.  This  makes  it  very  difficult  to 
establish  which  process  is  critical  to  the  consumption  of  NO  in  the  nozzle.  However,  it 
is  clear  that  the  rates  of  the  Zeidovich  reactions  must  play  a  critical  role,  as  well  as  the 
rates  of  oxygen  and  nitric  oxide  dissociation  and  recombination.  Although  N  atoms  play  a 
critical  role  in  the  Zeidovich  reaction  loop,  the  rate  of  N  atom  recombination  by  three-body 
processes  is  so  slow  that  it  does  not  play  a  role  in  this  flow.  Therefore,  we  consider  the 
sensitivity  of  the  nozzle  exit  conditions  to  these  critical  rates. 
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Axial  Distance  from  Throat  (m) 


Figure  18.  NO  formation/destruction  rate  due  to  recombination  and  Zeldovich  reactions 
on  nozzle  centerline. 

Recombination  of  Oxygen  Atoms 

A  critical  process  in  the  nozzle  flow  is  the  recombination  of  oxygen  atoms.  This  takes 
place  through  the  three-body  process  of  O  -f  O  -I-  M  — »  O2  +  M.  At  lower  temperatures, 
other  processes  may  play  a  role  in  this  process  (notably  O3).  As  discussed  above,  conven¬ 
tion  dictates  that  the  backward  reaction  rate  is  computed  from  the  forward  rate  via  the 
principle  of  detailed  balance.  In  equilibrium,  this  in  fact  must  hold,  but  there  is  no  rea¬ 
son  why  the  recombination  process  must  be  dictated  by  what  is  essentially  an  equilibrium 
thermodynamic  argument.  Rather,  the  recombination  process  can  proceed  through  many 
complicated  pathways  so  long  as  the  resulting  gas  composition  reaches  an  equilibrium 
condition. 

Keck  (1960)  and  Shui  et  al.  (1970)  develop  a  variational  theory  to  predict  three-body 
recombination  rates;  more  recently,  Korobeinikov  (1976)  developed  a  kinetic- theory  based 
approach  for  computing  three- body  recombination  rates.  These  theoretical  approaches 
compare  favorably  with  experimental  data.  Baulch  (1973,  1976)  compiled  the  available 
experimental  data  for  the  three- body  recombination  of  O  atoms  by  various  collision  part¬ 
ners.  Figure  19  shows  a  synthesis  of  their  results  for  the  relevant  collision  partners.  As  we 
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saw  from  the  previous  figures,  the  critical  recombination  process  occurs  when  the  gas  is 
between  about  3000  and  6000  K,  For  this  temperature  range,  there  is  about  a  factor  of  30 
variation  in  the  measured  recombination  rates;  in  general  the  most  active  collision  partner 
are  O  atoms,  and  the  least  active  N2  molecules.  The  theoretical  rates  fall  within  this 
general  scatter  of  the  experimental  data.  This  huge  variation  in  the  data  is  not  surprising 
given  the  difficulty  of  making  the  rate  measurements  and  the  various  systems  in  which 
they  were  made. 

Also  plotted  in  Figure  19  is  the  recombination  rate  based  on  detailed  balance  and 
the  dissociation  rates  used  in  our  work.  The  strange  variation  in  the  rate  at  temperatures 
less  than  1000  K  is  because  of  how  the  two  exponentially  varying  functions  for  the  forward 
reaction  rate  and  the  equilibrium  constant  interact  with  each  other  at  low  temperatures. 
This  behavior  is  not  reflected  in  the  data,  but  it  is  not  important  to  the  nozzle  flow  since 
the  flow  is  already  frozen  by  the  time  it  has  cooled  to  1000  K.  In  any  case,  it  is  clear  that 
our  recombination  rate  with  N2  as  collision  partners  falls  well  above  the  experimental  data 
at  the  critical  temperature  range.  The  rate  for  O-atoms  as  collision  partners  is  perhaps  a 
little  high,  but  within  the  general  range  of  the  data. 


Figure  19.  O  -f-  O  +  M  recombination  rates  from  Baulch  compilation  and  from  detailed 
balance  approach. 
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Therefore,  we  could  be  entitled  to  reduce  the  0  +  0  +  N2  — *  O2  +  N2  reaction  rate  by 
a  significant  margin  (say  by  a  factor  of  10).  However,  we  cannot  just  lower  the  rate  of  this 
process  because  then  the  reservoir  gas  would  no  longer  be  at  the  proper  thermodynamic 
equilibrium.  Thus,  we  must  lower  both  the  forward  and  backward  rate  of  this  reaction  to 
maintain  equilibrium  in  the  reservoir.  If  we  make  this  change  to  this  rate,  the  resulting 
nozzle  exit  plane  conditions  are  a  NO  mass  fraction  of  4.46%  (down  from  5.42%),  but  a 
lower  velocity  of  4317m/s  and  a  higher  O  atom  mass  fraction  of  3.10%.  Thus,  the  energy 
stored  in  the  chemical  energy  is  now  612kJ/kg,  up  by  a  factor  of  nearly  two.  Thus,  this 
change  moves  the  nitric  oxide  level  in  the  right  direction,  but  makes  the  velocity  worse. 

Zeldovich  Reaction  Rates 

Given  the  above  discussion  of  the  critical  role  of  the  Zeldovich  reactions  in  the  re¬ 
combination  of  air,  it  makes  sense  to  consider  the  effect  of  these  reaction  rates.  In  the 
present  work,  we  have  used  the  rates  of  Park:  kj  =  6.4  x  1017T_1  exp(-38400/T')  and 
kj  =  8.4  x  1012exp(-19450/T)  (cm3/moIos)  for  the  first  and  second  Zeldovich  reactions. 
However,  there  is  a  very  large  scatter  in  the  measured  rates.  (We  will  spare  the  reader 
additional  plots  showing  this  variation;  these  are  readily  available  in  the  Baulch  et  al. 
publications.)  Also,  quasi-classical  trajectory  calculations  show  that  the  Park  rates  are 
questionable  at  higher  temperatures  [Bose  et  al.  (1996,  1997)  These  calculations  using  ac¬ 
curate  potential  energy  surfaces  should  increase  in  accuracy  at  high  temperature,  contrary 
to  the  experimental  data.  Therefore,  let  us  consider  the  effect  of  using  these  QCT  rates  in 
the  nozzle  flow  field. 

We  have  run  several  cases  with  the  Bose  QCT  rates;  in  all  cases,  the  nitric  oxide  mass 
fraction  increases  at  the  nozzle  exit  plane.  Using  these  rates  for  both  reactions  gives  a 
NO  mass  fraction  of  6.70%;  this  is  primarily  due  to  the  slower  forward  rate  of  the  first 
Zeldovich  reaction.  Thus,  the  nozzle  flow  is  sensitive  to  these  rates,  it  does  not  appear 
that  accepted  rate  expressions  result  in  improved  agreement  with  the  Parker  data. 

Possible  Role  of  Excited  Electronic  States 

The  reservoir  temperature  (6318  K)  is  sufficiently  high  that  they  could  play  a  role  in 
the  recombination  process.  The  important  electronic  states  are  the  1 D2  state  of  O-atoms 
(2,0  eV  above  ground)  and  the  o1  and  bl  £+  states  of  O2  molecules  (1.0  and  1.6  eV 
above  ground).  All  of  these  states  are  metastables,  with  long  radiative  lifetimes.  There 
are  various  ways  in  which  these  species  could  play  a  role  in  the  recombination  process. 
However,  accepted  chemical  kinetics  mechanisms  and  associated  rates  are  not  available. 
Therefore,  we  must  postulate  possible  processes  and  assume  reasonable  rates  to  provide  a 
rough  bounds  for  their  relevance. 
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First  consider  first  excited  state  of  0  atoms;  denoting  this  state  as  O*,  it  is  possible 
that  the  second  Zeldovich  reaction  could  be  accelerated  by  this  species.  The  reaction 

NO  +  0*  — ►  O2  +  N 

is  now  exothermic  because  of  the  elevated  heat  of  formation  of  O’.  Presumably,  this 
reaction  would  be  very  fast,  resulting  in  the  rapid  destruction  of  NO  and  the  formation  of 
N  atoms  that  can  participate  in  the  first  Zeldovich  reaction.  However,  there  is  a  limited 
supply  of  O*  unless  there  is  a  means  of  generating  it  during  the  expansion.  The  only 
obvious  process  is 

0  +  0-f-0->02  +  0* 

which  would  yield  an  O*  atom  due  to  the  exothermic  recombination  of  02  by  O  atom  inter¬ 
actions.  This  reaction  has  the  potential  to  produce  a  significant  number  of  electronically 
excited  species. 

We  implemented  these  reactions  in  the  quasi-one-dimensional  code  and  reran  Park’s 
case.  To  do  so,  wc  must  assume  several  rate  constants,  and  assess  the  sensitivity  of 
the  results  to  the  choice  of  the  unknown  rate  constants.  Previously,  we  saw  that  the 
recombination  rate  of  02  via  ground  state  species  is  highly  uncertain;  clearly  the  rates 
that  we  need  for  these  reactions  are  even  more  uncertain.  Thus,  the  purpose  here  is  to 
see  if  this  process  has  any  potential  importance  in  the  nozzle  flow.  Here  we  take  the  most 
favorable  assumption  that  the  O*  species  do  not  radiatively  or  collisionally  decay;  rather 
they  are  only  destroyed  via  the  above  reactions.  Also,  to  provide  an  upper  bound,  we 
assume  that  all  0  +  O  +  O  reactions  result  in  the  formation  of  0*. 

In  order  to  implement  the  above  model,  we  need  an  expression  for  the  equilibrium 
constants  for  the  reactions.  We  use  statistical  mechanics  to  get: 

[O*]  _  9ie~B,^T 

[0]  ~  go+gie~°'/T 

with  £0  =  9,  pt  =  5,  and  0\  =  22830  K.  This  expression  can  then  be  used  to  adjust  the 
ground-state  equilibrium  constants  for  the  reactions  involving  O*. 

We  find  that  this  mechanism  does  not  appear  to  play  an  important  role  in  NO  de¬ 
struction.  There  is  only  a  small  number  (0.18%  mass  fraction)  of  O*  atoms  in  the  reservoir. 
Their  number  diminishes  quickly  as  the  gas  expands  due  to  the  very  strong  dependence  of 
the  equilibrium  excited  state  fraction  on  temperature.  Furthermore,  the  02  recombination 
process  does  not  form  many  excited  state  species,  again  because  of  the  extreme  tempera¬ 
ture  dependence  on  the  equilibrium  constant  for  this  reaction.  This  conclusion  appears  to 
be  independent  of  the  assumed  rates  for  these  reactions. 
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There  is  at  least  one  caveat  to  the  conclusion  that  the  0  excited  state  does  not  play  a 
role  in  the  destruction  of  NO.  We  have  assumed  a  simple  kinetics  model  for  0\  however 
if  there  were  additional  reaction  pathways,  it  may  be  possible  to  sustain  a  larger  level  of 
excited  atoms  within  the  nozzle.  However,  this  will  require  much  more  work. 

Now  let  us  turn  to  the  role  of  the  electronically  excited  states  of  O2.  During  the 
recombination  of  oxygen,  the  following  reactions  could  occur: 

O  +  O  +  M  -*■  Oa(X)  +  M 

0  +  0  +  M-»02(A)  +  M 
O  +  0  +  M  —  02{B)  +  M 

That  is,  the  oxygen  molecules  could  recombine  to  electronically  excited  states.  These 
species  would  store  additional  chemical  energy,  and  would  presumably  be  highly  reactive. 
Again,  we  face  the  issue  that  there  are  no  available  rate  data  for  these  reactions.  How¬ 
ever,  Keck  uses  his  variational  theory  approach  to  predict  these  rates;  he  shows  that  the 
recombination  to  the  A  state  is  about  a  factor  of  two  slower  than  to  the  ground  state,  and 
recombination  to  the  B  state  is  about  a  factor  of  ten  slower.  Thus,  there  does  not  appear 
to  be  any  favoring  of  these  reactions,  Furthermore,  these  reactions  would  store  additional 
chemical  energy  and  would  further  delay  the  relaxation  process.  As  such,  it  appears  that 
if  anything,  the  excited  02  states  would  tend  to  make  the  agreement  with  Parker’s  data 
even  worse. 

Additional  Comments 

Without  resorting  to  random  and  aphysical  modifications  of  the  chemical  kinetics 
model,  we  cannot  obtain  the  low  levels  of  NO  measured  by  Parker  in  the  LENS- 1  nozzle. 
There  could  be  additional  processes  that  we  have  not  considered.  Again,  it  should  be 
stressed  that  these  effects  must  be  quite  large. 

For  example,  there  is  N20  and  N02  in  the  reservoir  (at  levels  of  about  0.03%  to 
0.04%).  A  simple  kinetics  model  shows  that  these  species  are  rapidly  destroyed  during 
the  expansion,  and  do  not  play  a  role  in  the  destruction  of  NO.  We  have  neglected  the 
presence  of  ions  in  our  analysis;  NO+  is  the  most  prevalent  ion  at  8.1  x  10~5  mass  fraction. 
Tiiere  are  additional  trace  species  that  could  potentially  play  a  role,  though  again  these 
are  present  at  only  low  levels. 
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6.  Conclusions 

In  this  report  we  consider  the  flow  of  high-enthaipy  air  in  the  CUBRC  LENS-I  nozzle. 
We  compare  with  recent  noil-intrusive  laser  absorption  measurements  of  Parker  et  at. 
(2006,  2007)  in  an  attempt  to  better  understand  this  flow.  Parker’s  measurements  show  a 
high  velocity,  very  little  internal  energy  excitation,  and  a  low  (1.5%  by  mass)  level  of  nitric 
oxide  at  the  nozzle  exit  plane.  Our  calculations  predict  a  significantly  higher  level  of  NO 
(5.42%),  in  reasonable  agreement  with  more  detailed  simulations  of  Park  (2006).  We  also 
predict  a  correspondingly  low  velocity  and  a  higher  level  of  internal  energy  excitation  than 
Parker.  We  investigate  possible  reasons  for  this  large  discrepancy.  We  vary  the  Zeldovich 
rates  within  reasonable  bounds  from  the  literature,  and  we  use  experimental  data  to  modify 
the  O2  recombination  rate.  We  also  consider  the  possible  role  of  electronically  excited  O2 
and  O,  as  well  as  other  trace  species.  These  changes  make  some  differences  to  the  predicted 
NO  mass  fraction,  but  they  are  relatively  minor.  At  this  point,  we  cannot  reproduce 
Parker’s  measurements  using  kinetics  models  that  are  within  the  bounds  of  reality. 
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VI.  Summary 


During  the  past  three  years,  we  have  worked  to  develop  and  analyze  a  high-enthalpy 
flow  dataset  for  CFD  code  validation.  In  collaboration  with  CUBRC  Inc.  personnel,  we 
designed  the  experimental  conditions  for  the  double-cone  geometry  used  in  previous  low- 
enthalpy  nitrogen  tests.  We  then  used  our  standard  CFD  codes  and  thermo-chemical 
models  to  analyze  these  flows.  In  general,  the  comparisons  with  nitrogen  flows  are  good, 
even  at  high  enthalpy  (lOMJ/kg).  But  the  comparisons  between  predictions  and  experi¬ 
ments  for  air  above  5  MJ/kg  is  poor.  We  find  that  as  the  enthalpy  increases,  the  agreement 
gets  worse.  In  particular,  the  CFD  predicts  that  the  separation  zone  decreases  in  size  much 
more  rapidly  than  given  in  the  experiments. 

We  have  attempted  to  understand  the  source  of  this  discrepancy.  Based  on  our  pre¬ 
vious  work  with  nitrogen,  we  were  concerned  that  the  nozzle  flow  field  was  not  predicted 
correctly.  We  therefore  spent  a  lot  of  time  assessing  the  sensitivity  of  the  nozzle  test  section 
conditions  to  the  thermo-chemical  model  used  to  predict  the  nozzle  flow  field.  We  used 
a  vibrational  state-specific  model  to  assess  the  possibility  of  non- Boltzmann  vibrational 
energy  distributions  in  the  nozzle  -  this  was  completely  unsuccessful  in  explaining  the 
observed  discrepancies.  We  then  made  a  careful  study  of  the  assumptions  implicit  in  tiie 
thermo-chemical  model  for  air.  Again,  these  efforts  were  also  in  vain.  Therefore,  we  must 
conclude  that  there  is  another  unknown  explanation  for  the  discrepancy  between  the  mea¬ 
surements  and  the  simulations.  Future  work  will  involve  a  study  of  how  the  reservoir  flow 
behaves  during  the  shock  reflection  from  the  end-wall  of  the  shock  tunnel.  We  postulate 
that  the  reservoir  flow  does  not  reach  equilibrium  before  expanding  through  the  nozzle; 
clearly  this  would  change  the  facility  test-section  conditions. 
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